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Niels Bohr originally applied his approach to quantum mechanics to the H atom with great success. 
He then went on to show in 1913 how the same "planetary-orbit" model can predict binding for 
the H2 molecule. However, he misidentified the correct dissociation energy of his model at large 
internuclear separation, forcing him to give up on a "Bohr's model for molecules". Recently, we 
have found the correct dissociation limit of Bohr's model for H2 and obtained good potential energy 
curves at all internuclear separations. This work is a natural extension of Bohr's original paper 
and corresponds to the D — 00 limit of a dimensional scaling (D-scaling) analysis, as developed by 
Herschbach and coworkers. 

In a separate but synergetic approach to the two-electron problem, we summarize recent advances 
in constructing analytical models for describing the two-electron bond. The emphasis here is not 
maximally attainable numerical accuracy, but beyond textbook accuracy as informed by physical 
insights. We demonstrate how the interplay of the cusp condition, the asymptotic condition, the 
electron-correlation, configuration interaction, and the exact one electron two-center orbitals, can 
produce energy results approaching chemical accuracy. To this end, we provide a tutorial on using 
the Riccati form of the ground state wave function as a unified way of understanding the two- 
electron wave function and collect a detailed account of mathematical derivations on the exact 
one-electron two-center wave functions. Reviews of more traditional calculational approaches, such 
as Hartree-Fock, are also given. 

The inclusion of electron correlation via Hylleraas type functions is well known to be important, 
but difficult to implement for more than two electrons. The use of the D-scaled Bohr model off'ers 
the tantalizing possibility of obtaining electron correlation energy in a non-traditional way. 

PACS numbers: 31.10.-(-z, 31.15.-p, 31.25.-v, 31.50.-x 
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I. INTRODUCTION 



A. Overview 



We are in the midst of a revolution at the interface between chemistry and physics, largely due to the interplay 
between quantum optics and quantum chemistry. For example, the explicit control of molecules afforded by modern 
femtosecond lasers and adaptive computer feedback 1] has opened new frontiers in molecular science. In such studies, 
molecules are controlled by sculpting the amplitude and phase of femtosecond pulses, forcing the molecule into 
predetermined electronic and rotational-vibrational states. This holds great promise for vital applications, from the 
trace detection of molecular impurities, such as dipicolinic acid as it appears in anthrax 0, to the utilization of 
molecular excited states for quantum information storage and retrieval [2|. 

We are thus motivated to rethink certain aspects of molecular physics and quantum chemistry especially with regard 
to the excited state dynamics and coherent processes of molecules. The usual discussions of molecular structure are 
based on solving the many-particle Schrodinger equation with varying degree of sophistication, from exacting Diffusion 
Monte Carlo methods, coupled cluster expansion, configuration interactions, to density functional theory. All are 
intensely numerical. Despite these successful tools of modern computational chemistry, there remains the need for 
understanding electron correlations in some relatively simple way so that we may describe excited states dynamics 
with reasonable accuracy. 

In this work, we propose to reexamine these questions in two complementary ways. One approach is based on the 
recently resurrected Bohr's model for molecules In particular, we show that by modifying the original Bohr's 

model [3 in a simple way, specially when augmented by dimensional scaling (D-scaling), we can describe both the 
singlet and triplet potential of H2 with remarkable accuracy (see Figs. |5J|5l. 

In another approach, following the lead of the French school of Le Sech we use correlated two-center orbitals 

of the molecule to model H2's ground and excited state. This approach worked well, even when only a simple 
electron correlation function is used, see Table 



TABLE I: The binding energy of H2 molecule based on "exact" two-center orbitals. 



Orbital 


Binding energy (eV) 






No Free Parameter 


1 Free (screening) 






parameter 


Jatte lYI-WJ 


4.50 


4.60 


HvUeraas IVI.SSl 4.51 


4.62 



When we allow the a and B parameters of Eqs. ljVI.55|l and ljVI.58|l to vary, we obtain a binding energy of 4.7 eV. The 
binding energy is comparable to the experimental value of 4.7 eV. 

The Bohr model and D-scaling technique taken together with good (uncorrelated) molecular orbitals is especially 
interesting and promising. As discussed in subsection C, the Bohr model yields a good approximation to the electron- 
electron Coulomb energy, which can be used to choose a renormalized nuclear charge and a much improved (correlated) 
two electron wave function. 



B. The Bohr molecule 



Figure n displays the Bohr model for a hydrogen molecule ^ IE l^i in which two nuclei with charges Z\e\ are 
separated by a fixed distance R (adiabatic approximation) and the two electrons move in the space between them. 
The model assumes that the electrons move with constant speed on circular trajectories of radii pi = p2 = P- The 
circle centers lie on the molecule axis z at the coordinates zi = ±Z2 = z. The separation between the electrons is 
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FIG. 1: Cylindrical coordinates (top) and electronic distances (bottom) in H2 molecule. The nuclei Z are fixed at a distance K 
apart. In the Bohr model, the two electrons rotate about the internuclear axis z with coordinates p\, z\ and p2, 22 respectively; 
the dihedral angle between the (pi, z\) and (p2, ^2) planes remains constant at either <^ = tt or <^ = 0. The sketch corresponds 
to configuration 2 of Fig. |51 with <j!> = tt. 




constant. The net force on each electron consists of three contributions: attractive interaction between an electron 
and the two nuclei, the Coulomb repulsion between electrons, and the centrifugal force on the electron. We proceed 
by writing the energy function _B = r + where the kinetic energy T = p^/2m + p|/2m for electrons 1 and 2 can be 
obtained from the quantization condition that the circumference is equal to the integer number n of the electron de 
Broglie wavelengths 27rp = nh/p, so that we have T — /2m — v?T? jlmp^; the unit of distance is taken to be the 
Bohr radius oq = /me^, and the unit of energy the atomic energy, e^/ag where m and e are, respectively, the mass 
and charge of the electron. The Coulomb potential energy V is given by 

Tal ni ra2 n2 ri2 R 

where Tai {i = 1,2) and r^i are the distances of the ith electron from nuclei A and B, as shown in Fig. ^ (bottom), 
ri2 is the separation between electrons. In cylindrical coordinates the distances are 

ni = 






here R is the internuclear spacing and (j) is the dihedral angle between the planes containing the electrons and the 
internuclear axis. The Bohr model energy for a homonuclear molecule having charge Z is then given by 

^^=U4 + 4)+nPl,P2,Zl,Z2,0). (1.2) 

^ \ Pi P2 J 

Possible electron configurations correspond to extrema of the energy function (|I.2|I . For ni = 712 = 1 the energy has 
extrema at pi — P2 = P, zi = ±Z2 = z and = tt, 0. These four configurations arc pictured in Fig. |2] (upper panel). 
For example, for configuration 2, with zi = —Z2 = z, 4> = ir, the extremum equations dE/dz = and dE/dp = read 

Z{R/2-z) z 

[p2 + (R/2 - z)2]3/2 4[^2 + ^2]3/2 



Z{R/2 + z) 
[p2 + (i?/2 + z)- 



13/2 



(1.3) 
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Zp 



Zp 



4[p2 + ;22]3/2 p3 



(1.4) 



which are seen to be equivalent to Newton's second law applied to the motion of each electron. Eq. (|I.3p specifies that 
the total Coulomb force on the electron along the z— axis is equal to zero; Eq. HI.4|I specifies that the projection of 
the Coulomb force toward the molecular axis equals the centrifugal force. At any fixed internuclear distance i?, these 
equations determine the constant values of p and z that describe the electron trajectories. Similar force equations 
pertain for the other extremum configurations. 

In Fig. 121 (lower panel) we plot E(K) for the four Bohr model configurations (solid curves), together with "exact" 
results (dots) obtained from extensive variational calculations for the singlet ground state and the lowest triplet 
state, '^Ej m. In the model, the three configurations 1, 2, 3 with the electrons on opposite sides of the internuclear 
axis ((/) = tt) are seen to correspond to the singlet ground states, whereas the other solution 4 with the electrons 
on the same side (0 = 0) corresponds to the first excited, triplet state. At small internuclear distances, the 
symmetric configuration 1 originally considered by Bohr agrees well with the "exact" ground state quantum energy; 
at larger _R, however, this configuration's energy rises far above that of the ground state and ultimately dissociates 
to the doubly ionized limit, 2H++2e. In contrast, the solution for the asymmetric configuration 2 appears only for 
i? > 1.20 and in the large R limit dissociates to two H atoms. The solution for asymmetric configuration 3 exists only 
for R > 1.68 and climbs steeply to dissociate to an ion pair, H++H^. The asymmetric solution 4 exists for all R and 
corresponds throughout to repulsive interaction of two H atoms. 



1 symmetric (Bohr) 
Zi=Z2=0, (|)=Jl 



2 asymmetric 
Zi=-Z2 , (|)=7l 



3 asymmetric 

Zi=Z2 , <^—n 



4 asymmetric (repulsive) 
Zi=-Z2 , <|)=0 



-0.4 
-0.5 
-0.6 
-0.7 
-0.8 
-0.9 
-1.0 
-1.1 ■ 
-1.2- 




H^^H^-i-H" (R>1.68) 



• 2 H2^2H (R>1.2) 
"exact" values 



3 4 

R, a.u. 
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FIG. 2: Energy E{R) of H2 molecule for four electron configurations (top) as a function of internuclear distance R calculated 
within the Bohr model (solid lines) and the "exact" ground and first excited '^Ej state energy of Ref. Q (dots). Unit of 
energy is 1 a.u.= 27.21 eV, and unit of distance is the Bohr radius. 

We then extend these "Bohr molecule" studies in several ways. In particular, we use a variant of the dimensional 
scaling (D scaling) theory as it was originally developed in quantum chromodynamics and applied with great success 
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to molecular and statistical physics |10|, [ll] . This is based on an analysis in which the usual kinetic energy terms in 
the Schrodinger equation are written in D dimensions, i.e., 



(15) 

2m ^ dxf ^ 2m ^ dx^ ' 

This provides another avenue into the interface between the old (Bohr-Sommerfeld) and the new (Heisenberg- 
Schrodinger) quantum mechanics. In particular, when D oo the two electron Schrodinger equation can be scaled 
and sculpted into a form which is when Z) — > cxd identical to the Bohr theory of the H2 molecule, and much better 
than when l/D or other corrections are included. 

C. Simple correlation energy from the Bohr model 

The Bohr model offers an effective way to treat most of the correlation energy absent in the conventional Hartree- 
Fock (HF) approximation. Here we show how a charge renormalization method can be applied to improve the ground 
state energy obtained in the HF approximation. We start from He-likc ions and consider a nucleus with charge Z and 
two electrons moving around it. According to the Bohr model the ground state energy is given by the minimum of 
the following expression 

1 f 1 1\ Z Z 1 

E ^ 7: { ^ + ^] + , (1-6) 

2 \Pi Pi J Pi P2 ^JpI+pI - 2pip2 COS0 



where (f> is the dihedral angle between electrons. At the minimum cj) = n and Eq. Ijl.6|l reduces to 



1 f I 1\ Z Z I , , 

E = -\ ^ + ^] + ■ (1-7) 

2 \Pi pi J Pi P2 P1+P2 



Optimization with respect to pi and p2 yields 



^.2 



P^^P^ = J-T^ EB = -[Z--j . (1.8) 

The HF approximation in the framework of the Bohr model means that optimum parameters pi and p2 are determined 
by minimization of Eq. 11. 7() with no electron repulsion term, i.e., omitting the electron-electron correlation. In the 
HF approximation the Bohr model gives 

Pl=P2 = ^, ^B-HF = -^' + |. (1.9) 

For the He atom {Z = 2) we obtain E-q-uf — —3 a.u., while Eb — —3.0625 a.u. Thus, the inclusion of correlation 
shifts the ground state energy down by 

E^^^^ " " -0.0625 a.u. (1. 10) 

The Bohr model itself is quasiclassical and, as a consequence, it predicts the He ground state energy with only 5.4% 
accuracy (-Eoxact = —2.9037 a.u.). However, the Bohr model provides a quantitative way to include the correlation 
energy. Let us consider the He ground state energy calculated using the HF (effective charge) variational wave function 

*(ri,r2)=Cexp[-Z(ri+r2)], (1.11) 

where Z is a variational parameter (effective charge), which is determined by minimizing the energy E ^ Z^ — 2ZZ -\- 
5Z/8, Z = Z — 5/16. The wave mechanical HF energy is 



<F = -(^-^) . (1.12) 
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TABLE II: Ground state energy of the He-like ions in the HF approximation -Ehf ^^'^ the value improved by the Bohr model 
Eb- The last two columns compare the accuracy of the HF and the improved result. 



7 

Zj 




7 7 a 




J-^cxact 






9 






9 Qi ni 


9 00*^7 


1 0"^ 


n 99 


Q 
O 




u.uouo 




7 97QQ 


n 7Q 


n 079 


4 


-13.597 


0.0540 


-1 3 6609 


-1 3 6555 


0.42 


0.033 


5 


-21.9726 


0.0558 


-22.0351 


-22.0309 


0.26 


0.019 


6 


-32.3476 


0.0570 


-32.4102 


-32.4062 


0.18 


0.012 


7 


-44.7226 


0.0578 


-44.7852 


-44.7814 


0.13 


0.008 


8 


-59.0976 


0.0584 


-59.1602 


-59.1566 


0.10 


0.006 


9 


-75.4726 


0.0589 


-75.5352 


-75.5317 


0.08 


0.004 


10 


-93.8476 


0.0592 


-93.9102 


-93.9068 


0.06 


0.003 



For Z = 2 we obtain E^p — —2.8476 a.u. The difference between E^p and the exact value is due to the correlation 



energy missing in the HF treatment. One can notice that if we add the correlation energy Ijl.lOp to E^p we obtain 

EnF + Eforr - -2.9101 a.u.. 



which substantially improves the answer and deviates by only 0.2% from the exact value. Such an idea can be 
incorporated by renormalization of the nuclear charge [l^ . Let us define an effective charge Z^ff by the condition 

£^B-HF(^cff) = E^p{Z), (1-13) 

which yields 




(1.14) 

The effective charge improves the Bohr model energy by taking into account the difference between the quasiclassical 
and fully quantum mechanical description. The effective charge is calculated from the correspondence between the 
Bohr model in the HF approximation and the quantum mechanical HF answer. Now, if we take the Bohr model 
energy Eb{Z) 11. 8|) (that includes correlation) but with Z^s instead of Z it improves the quantum mechanical HF 
answer: 



Es{Z.,,) = -[Z-^\ -1- = EZ{Z)-y^. (1.15) 



The correction energy —1/16 is independent of Z and coincides with Eq. Ijl.l0|) . Table iHl compares the quantum 
mechanical HF answer for the ground state energy of He- like ions and the improved value (|I.15|) . Depending on Z 
Eq. (|I.15|I improves the accuracy 10 — 20 times. 



D. Correlated two-center orbitals 



^From the preceding discussion it is clear that we need good (hopefully simple) HF wave functions. There is, of 
course, a g reat deal of work on this problem but we find the two-center orbital approach of Le Sech and coworkers QilS 
IT^ IT3 . Hsf and of Patil 16] to be especially useful. In a previous publication 17], we attempted a first principle (semi- 
tutorial) presentation employing that the exact two-center orbitals obtained from solving the Schrodinger equation 
for the HJ ion. As shown by Le Sech these are the most useful building blocks for constructing the electronic wave 
functions of the homonuclear H2 molecule. One simple form of the electronic ground state constructed with two-center 
orbitals is 

vE'H.(l,2) = 4'H+,ia(l)*H+,ia(2)Xoo {l + \ri2^ , (1.16) 

where ^1^^+ is the solution of the Schrodinger equation for the HJ in prolate-spheroidal (ellipsoidal) coordinates, 
Xoo = [| Tii2) — 1 iit2)]/\/2 is singlet spin function, and [l + jri2) is the Hylleraas correlation factor. See more detailed 
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discussions of (|I.16|) in Sections IV- VI below. This wave function yields a binding energy of 4.5 eV for H2 molecule 
without any variational parameters. Variation with respect to a couple of parameters in the function HI.16|I shifts 
the binding energy to 4.7 eV, giving remarkable agreement with the experimental value. To achieve the same result, 
sums over many one-centered atomic orbitals or Hylleraas type of wave function (cf. (|IV.30|I below) that explicitly 
include the interelectronic distance are usually used. This has been demonstrated by the earlier work of Kolos and 
Wolniewicz 0|. Kolos and Szalewicz and that of James and Coolidge 0|, respectively. 

In these studies the introduction of a correlation factor, taking into account the Coulomb interaction between the 
two electrons, is naturally motivated by considering the trial wave function as broken into three parts, we write 

«'(ri,r2) = *(ri)^'(r2)/(ri,r2), (1.17) 

where 4'(ri) and 5'(r2) are exact one electron solutions in the absence of interaction between electrons. For ^'(ri, r2) 
the Schrodinger equation (in atomic units) gives 

*(r2)/(ri,r2) f-^ _ ^ - ^) ^(r^) + *(ri)/(ri, r2) ( - ^ - ^) $(.2) 
V 2 ral ni J V 2 ra2 n2 J 

+ ^'(ri)^'(r2) — -i - ^ H f{ri,r2) + cross terms 

V 2 2 ri2 J 

E-^^^{n)^{r2)fin,r2). (1.18) 

where cross terms mean terms that go as Vi^'(ri) • Vi/(ri, r2), etc. The solution with only the first two terms is just 
that for H^. The functions ^'(ri) and ^'(r2) exponentially decay at large distances from the nuclei. The third term 
corresponds to the Schrodinger equation for two free electrons, 

V2 V2 1 \ 

/(n,r2) = £/(ri,r2) (1.19) 



2 2 ri2 

that is the solution to Eq. (|I.19|) is well known and is given in terms of Coulomb wave functions, i.e., confluent 
hypergeometric functions (see more detailed discussions in Subsection V.B), which yields the (1 + fi2/2) Hylleraas 
factor as an asymptotic at small ri2. In order to place this part of the present review in perspective and be ready for 
the paradigm shift from the old quasi quantum-mechanical Bohr model to the new fully wave-mechanical Schrodinger- 
Born-Oppcnhcimcr model, we next give a brief history of the molecular orbital concepts and computations. 



E. Context 



Molecular quantum chemistry is a fascinating success story in the annals of 20th Century science. In 1926, 
Schrodinger introduced the all-important wave equation which soon bore his name. In the following year Schrodinger's 
new theory was applied to the simplest molecular systems of the hydrogen molecular ion by Burrau [l^l and to 
the hydrogen molecule H2 by Heitler and London [2J| and Condon |2!tI |. In the same year. Born and Oppenheimer l2a| 
published their important paper dealing with molecular nuclear motion. Further, in 1928, Hund and Mulliken |27| 
presented their venerable molecular orbital (MO) theory, which provided a powerful computational tool for chemistry 
and a foundation for the subsequent development of modern molecular science. 

Diatomic molecules such as H2 and HeII+ are the simplest of all molecules. Their analysis, modeling and compu- 
tation constitute the bedrock of the study of chemical bonds in molecular structures. To quote a recent insightful 
article by Cotton and Nocera 

"It can be said without fear of contradiction that the two-electron bond is the single most important 
stereoelectronic feature of chemistry." 

Indeed, the description of the covalent bond in diatomics, based on the methods of Heitler-London and Hund-MuUiken, 
is one of the crowning achievements of quantum mechanics and fundamental physics. 

Computational quantum chemistry dawned in 1927 with the advent of the Heitler-London method. However, the 
accuracy of these early numerical results were not satisfactory, as can be seen from the following quotation (Hinchliffe 
m p. 254]): 

"The calculated bond length and dissociation energy [of MO theory] are in poor agreement with experiment 
than those obtained from simple VB [valence bond] treatment (Table 15.3), and this puzzled many people 
at the time. It has also led them to believe that the VB method was the correct way forward for the 
description of molecular structure ... ." 
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New ideas were then proposed to improve the numerical accuracy of the Heitler-London and Hund-Mulhken method. 
The first idea of configuration interaction (CI) is to incorporate excited states into the wave function. The second 
idea of correlation introduces exphcit dependence of the interelectronic distance in the wave function. The idea of 
correlation was first demonstrated by Hylleraas 's^ in 1929 for the helium atom and by James and Coolidgc in 
1933 for H2. The use of configuration interaction and correlation are key evolutionary steps in improving the original 
ideas of Heitler-London and Hund-MuUiken. We can quote the following from Rychlewski '■ 

"... Very soon it has been realized that inclusion of interelectronic distance into the wave function is 
a powerful way of improving the accuracy of calculated results .... Today, methods based on explic- 
itly correlated wave functions are capable of yielding the "spectroscopic" accuracy in molecular energy 
calculations (errors less than the orders of one /i Hartree) ... ." 

For more than two electrons, it is difficult for most numerical methods to include electron correlations directly except 
in Monte Carlo simulations. When it is possible, as in the two electron case, excellent results can be achieved with 
very compact wave functions. 

Molecular calculations are inherently more difficult than atomic calculations. The fundamental difficulty is well 
stated by Teller and Sahlin [3^: 

"The molecular problem has a greater inherent complexity than the corresponding atomic problem.... In 
atoms, degeneracy due to spherical symmetry causes many levels to have nearly the same energy. This 
grouping of levels is responsible for the presence of a shell structure in atoms, and this shell structure is 
in turn the primary reason for the striking and simple behavior of atoms and the consequent successes 
of the independent-electron approximation for atomic systems. In passing from atoms to molecules the 
symmetry is reduced and the amount of degeneracy for the electronic levels becomes smaller, and, as a 
consequence, the power of the independent-particle model is decreased relative to the atomic case." 

Nothing illustrates this loss of symmetry, and its consequent loss of validity of the independent-particle picture better 
than the complete failure of the molecular orbital picture to account for the correct dissociation energy of II2 . At large 
internuclear separation, the symmetry is greatly reduced, and the independent occupation of single-particle molecular 
orbitals fails catastrophically. This failure can of course be averted by configuration-interaction, but this extra work 
makes it plain that molecular problems are inherently more complicated than atomic problems. Fortunately, for the 
investigation of ground and excited molecular states near equilibrium, one is far from the dissociation limit; the loss 
of symmetry complicates the calculation of the molecular orbital, but the independent particle model remains a good 
approximation . 

In the case of II2 , a natural candidate is the orbital of the two- center one- electron molecular ion. Such orbitals will be 
referred to as diatomic orbitals (DO) or, in more complicated cases, shielded diatomic orbit als ( SDO) when shielding is 
a factor to be considered. The early (1930s) ansatz wave functions of James and Coolidge are expressed in terms 
of prolate spheroidal coordinates of the two electrons with respect to the two centers of the diatomic nuclei. However, 
their wave functions are not DOs in that they are not expansions of the exact one electron Hj states. Rather, their 
approach is CI with a basis conveniently chosen for numerical evaluation. Their work is the forerunner of the Polish 
quantum chemistry group [ol llSlIT^ of Kolos, Wolniewicz, etc., which have achieved the highest accuracy in numerical 
computation of two-electron molecules. The high accuracy obtained by Kolos and Wolniewicz in is admirable, 
but as noted by Patil, Tang and Toennies p^ . 

". . . It is, however, perhaps somewhat unfortunate that these very impressive accomplishments have largely 
discouraged further fundamental studies on novel approaches to obtain accurate wave functions more 
directly ... ." 

A similar comment was made much earlier by Mulliken |33j |: 

"[T]he more accurate the calculations become the more concepts tend to vanish into thin air." 

Thus the human quest for comprehension remains, and the recent research on novel approaches to obtain accurate 
wave functions have indeed yielded accurate, physically motivated, and compact two-electron wave functions. Patil 
et al. |34l l35l l37l l38l l39l |40| have advocated the construction of coalescense wave functions by incorporating both 
cusp and asymptotic conditions. We have provided a detailed review with simplfied derivations of this development 
in Section HI. The other approach is the use of diatomic orbitals. Historically, the original idea of using DOs as 
basis for diatomic molecules seems to begin from the work of Wallis and Hulburt ^| . More extensive references and 
history can be found in the works of Mclean, Weiss and Yoshimine 42], Teller and Sahlin [s^l and ShuU ^3[. Wallis 
and Hulburt's result was not particularly successful, because there was no explicit electron correlation and solving 
the two-center wave function was difficult. According to Aubert, Bessis and Bessis 0, Part I, p. 51]: 
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". . . the use of these functions, i.e., diatomic orbitals (DOs), within the one-configuration molecular-orbital 
scheme has not been very successful, owing to the difficulty of taking into account the interelectronic 
interactions and, moreover, owing to the complexity of calculations." 

The calculation of wave function improved over the years, culminating in the extensive tabulations by Teller 
and Sahlin 32]. In 1974-75, Aubert, Bessis and Bessis published a three-part series 13] detailing how to determine 
SDOs for diatomic molecules. These three papers emphasize the determination of shielded DOs. Sur pris ingly, the use 
of DO with correlation to study H2 was not undertaken until 1981 by Aubert-Frecon and Le Sech [Tj]. Le Sech, et 
al. have since then made further refinements to the method (Siebbeles and Le Sech Le Sech 0). 

Our study of the DO's approach was motivated by our strong interest in the modeling and computation of molecules. 
We were especially attracted by DOs as a natural and accurate description of chemical bonds. In Scully et al. [T^ . 
largely unaware of the prior work done by Aubert-Frecon and Le Sech of the French school, we obtained simple 
correlated DOs for diatomic molecules with good accuracy. The present paper represents part of our continued 
efforts in this direction. In this work, we will first study the mathematical properties of wave functions such as their 
cusp conditions, asymptotic behaviors, and forms of correlation functions. We summarize methods, techniques and 
formulas in a tutorial style, interspersed with some unpublished results of our own. It it not our intention to complete 
an exhaustive review on this vast subject, rather, only to record developments relevant to our interest in sufficient 
details. We apologize in advance for any inadvertent omissions. 

F. Outline 

The present paper is organized as follows: 

(i) In Section ^ we present some recent progress of an interpolated Bohr model. 

(ii) In Section [nil we discuss some general and fundamental properties of atoms and molecules, including the Born- 
Oppenheimer separation, the Feynman-Hellman Theorem, Riccati form of the ground state wave function, 
proximal and asymptotic conditions, the coalescent construction and the dissociation limit. 

(iii) In Section IIVI wc introduce the basics of the 1-electron two-centered orbitals from the classic explicit solutions 
of Hylleraas and Jaffc. The one-centered and multi-centered orbitals will also be reviewed. 

(iv) In Section we present the details of the derivations of the all-important cusp conditions of Kato, and show 
examples as to how to verify them in prolate spheroidal coordinates. We also provide a glossary of various 
correlation functions that satisfy the interelectronic cusp condition. 

(v) In Section lvTI we discuss numerical modeling of diatomic molecules and compare results with the classic methods 
such as the Heitler-London, Hund-Muliken, Hartree-Fock and James-Coolidge, and the new approach of two- 
centered orbitals by Le Sech, et al. 

(vi) In Section FVIII we discuss alternative methods for molecular calculations based on the generalized Bohr model 
and the dimensional scaling. 

(vii) Finally, in Section IVlIII we give our conclusions and present an outlook. 

II. RECENT PROGRESS BASED ON BOHR'S MODEL 

The diatomic molecules in a fully quantum mechanical treatment addressed in Subsection LE requires solution of the 
many-particle Schrodinger equation. However, such an approach also requires complicated numerical algorithms. As a 
consequence, for a few electron problem the results become less accurate and sometimes unreliable. This is pronounced 
for excited electron states when the application of the variational principle is much less involved. Therefore, invention 
of simple and, at the same time, relatively accurate methods of molecule description is quite desirable. 

In this section we discuss a method which is based on the Bohr model and its modification |^ . In particular, 
we show that for H2 a simple extension of the original Bohr model @ describes the potential energy curves for the 
lowest singlet and triplet states just about as nicely as those from the wave mechanical treatment. 

The simplistic Bohr model provides surprisingly accurate energies for the ground singlet state at large and small 
internuclear distances and for the triplet state over the full range of R. Also, the model predicts the ground state 
is bound with an equilibrium separation i?e ~ 1-10 and gives the binding energy as Eb ~ 0.100 a.u.= 2.73 eV. The 
Heitler-London calculation, obtained from a two-term variational function, yields Re = 1.51 and Eb = 3.14 eV [23 ]. 
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whereas the "exact" results are i?e = 1-401 and Eb = 4.74 eV 0- For the triplet state, as seen in Fig. [S] the 
Bohr model gives a remarkably close agreement with the "exact" potential curve and is in fact much better than the 
Heitler-London result (which, e.g., is 30% higher for i? = 2). One should mention that in 1913, Bohr found only the 
symmetric configuration solution, which fails drastically to describe the ground state dissociation limit. 

The simple Bohr model offers valuable insights for the description of other diatomic molecules. For N electrons the 
model reduces to finding configurations that deliver extrema of the energy function 

1 ^ 

^= oE^ + ^('"i''"2'-'r^'^)' (III) 

where the first term is the electron kinetic energy, while V is the Coulomb potential energy. In such formulation 
of the model there is no need to specify electron trajectories nor to incorporate nonstationary electron motion. Eq. 
(jlLip assumes that only at some moment in time the electron angular momentum equals an integer number of h and 
the energy is minimized under this constraint. In the general case, the angular momentum of each electron changes 
with time; nevertheless, the total energy remains a conserved quantity. 

Let us now discuss the ground state potential curve of HeH. To incorporate the Pauli exclusion principle one can 
use a prescription based on the sequential filling of the electron levels. In the case of HeH the three electrons can 
not occupy the same lowest level of HeH++. As a result, we must disregard the lowest potential energy curve E{R) 
obtained by the minimization of Eq. (|II.ip and take the next possible electron configuration, which is shown in Fig. 
01 (upper panel). For this configuration, ni ~ n2 = = 1, however, the right energy corresponds to a saddle point 
rather than to a global minimum. In order to obtain the correct dissociation limit we assign the helium an effective 
charge — 1.954. The charge matches the difference between the exact ground state energy of the He atom and 
the Bohr model result. Fig. |31 shows the ground state potential curve of HeH in the Bohr model (sohd curve) and the 
"exact" result (dots) obtained from extensive variational calculations 0|. The Bohr model gives a remarkably close 
agreement with the "exact" potential curve. 
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FIG. 3: Energy E{R) of HeH molecule for the electron configuration shown on the upper panel as a function of internuclear 
spacing R calculated within the Bohr model for ni — n2 = n-^ — 1, Z]^^ — 1.954 (solid line) and the "exact" ground state 
energy of Ref . ,4^ (dots) . 
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A. Interpolated Bohr model 

The original Bohr model assumes quantization of the electron angular momentum relative to the molecular axis. 
This yields a quite accurate description of the H2 ground state E{R) at small R. However, E{R) becomes less accurate 
at larger internuclear separation as seen in Fig. |21 To obtain a better result one can use the following observation. 
At large R each electron in H2 feels only the nearest nuclear charge because the remaining charges form a neutral H 
atom. Therefore, at large R the momentum quantization relative to the nearest nuclei, rather than to the molecular 
axis, must yield a better answer. This leads to the following expression for the energy of the H2 molecule 

E^oi^ + ^j + — + 17 ("-2) 

2 V^al K2J ''al ^'bl ra2 n2 ri2 R 

For ni = n2 = 1 and R > 2.8 the expression (III.2p has a local minimum for the asymmetric configuration 2 of Fig. 
121 We plot the corresponding E{R) without the 1/R term in Fig. ^ (curve 2). At i? < 2.8 the electrons collapse into 
nuclei. At small R we apply the quantization condition relative to the molecular axis which yields curve 1 in Fig. 0] 
To find E{R) at intermediate separation we connect smoothly the two regions by a third order polynomial (thin line). 
Addition of the 1/R term yields the final potential curve, plotted in Fig. [51 The simple interpolated Bohr model 
provides a remarkably close agreement with the "exact" potential curve over the full range of R. 
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FIG. 4: The Bohr model E{R) for H2 molecule without 1/R term. Curves 1 and 2 are obtained based on the quantization 
relative to the molecular axis (small R) and the nearest nuclei (large R) respectively. Thin line is the interpolation between 
two regions. 

As an example of application of the interpolated Bohr model to other diatomic molecules, we discuss the ground 
state potential curve of LiH. The Li atom contains three electrons, two of which fill the inner shell. Only the outer 
electron with the principal quantum number n = 2 is important in forming the molecular bond. This makes the 
description of LiH similar to the excited state of H2 in which two electrons possess ni = 1 and 712 = 2. So, we start 
from the H2 excited state and apply the interpolated Bohr model as described above. Then, to obtain E{R) for LiH, 
we take the H2 potential curve and add the difference between the ground state energy of Li (-7.4780 a.u.) and H in 
the n = 2 state, i.e., add —7.3530 a.u. The final result is shown in Fig. (lower solid line), while dots are the "exact" 
numerical answer from ;4:5j. One can see that the simple interpolated Bohr model provides quite good quantitative 
description of the potential curve of LiH, which is already a relatively complex four electron system. 

III. GENERAL RESULTS AND FUNDAMENTAL PROPERTIES OF WAVE FUNCTIONS 



Having examined the quasi-quantum mechanical Bohr model in the preceding two sections, we now attend to the 
fully quantum mechanical model and its approximation and analysis. 
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R, a.u. 

FIG. 5: Ground state E{R) of H2 molecule as a function of internuclear distance R calculated within the interpolated Bohr 
model (solid line) and the "exact" energy of Ref. Q (dots). 



A. The Born—Oppenheimer separation 



The underlying theoretical basis for problems involving a few particles in atomic and molecular physics is the 
Schrodinger equation for the electrons and nuclei, which provides satisfactory explanations of the chemical, electro- 
magnetic and spectroscopic properties of the atoms and molecules. Assume that the system under consideration has 

Ni 

Ni nuclei with masses AIk and charges cZk, e being the electron charge, for i^T = 1, 2, . . . , iVi, and N2 = J2 is 

K=\ 

the number of electrons. The position vector of the Xth nucleus will be denoted as while that of the /cth electron 
will be Tfe, for fc = 1, 2, . . . , iV2. Let m be the mass of the electron. The Schrodinger equation for the overall system 
is given by 



Ni 

E 

K=l 



2Mk 



K 



N-2 

E 

k=l 



2m 



Ni N2 

EE 

K=l k=l 



Rk -rk\ 



1 



N2 

E 



k^k' ' 



1 ^ Zk 



K,K' = 1 



R 



K' 



-^{R^r) = E'^iR,r), 



(III.1) 



where H is the Hamiltonian and 



R= {Ri,R2,. ■ ■ ,Rni), r = (ri,r2, . . . ,rArJ. 



The above equation is often too complex for practical purposes of studying molecular problems. Born and Oppen- 
heimer |26j | provide a reduced order model by approximation, permitting a particularly accurate decoupling of the 
motions of the electrons and the nuclei. The main idea is to assume that ^ in ljIII.l|) takes the form of a product 



«'(ii,r-) = G{R)F{R,r). 



(111.2) 
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FIG. 6; Ground (^S''') state energy E{R) of LiH molecule as a function of internuclear distance R calculated within the 
interpolated Bohr model (lower solid line) and the "exact" energy of [4fi| (dots). Upper solid curve is the first excited (^S^) 
state energy of LiH obtained from the Bohr model E{R) of H2 molecule by adding the difference between the ground state 
energy of Li and H. 



Substituting l|TTr2|) into (|ITlT1) . we obtain 
GiR) 



2 ^2 N2 



k=l 



1 



N2 

E 



k,k' = l 



F{R,r) 



where 



E 

K=l 



2Mk 



-V 



K 



1 



ZkZk'g" 



2 j^K' \^k-Rk' 

K,K' = 1 



G{R))F{R,r) 



%+% = EG{R)F{R,r), 



Ti ^ ~G{R) J2 ^^kG{R) ■ ^KF{R,r), 
T.^-G{R)Y^-—^lF{R,r). 



K=l 



2Mk 



(III.3) 



(III.4) 



(III.5) 



The essential step in the Born-Oppcnheimer separation consists in dropping Ti and T2 in ljIII.3|) . This leads to the 
separation of the electronic wave function 



/ 



-2 ^2 



N2 



2m 



Ev^n E 



Ni N2 

^rrr ~ E E 



k=l 



Xk — r, 1 

k^k' ' *- x=ife=i 

k,k' = l 



Rx-Vkl 



F{R,r) ^ E,{R)FiR,r), 



(III.6) 
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and a second equation for the nuclear wave function 



G{R) = EG{R), 



where Ee{R) is the constant of separation. 

In "typical" molecules, the time scale for the valence electrons to orbit about the nuclei is about once every 10~^^ 
s (and that of the inner-shell electrons is even smaller), that of the molecular vibration is about once every 10~^^ s, 
and that of the molecule rotation is every 10~^^ s. This difference of time scale is what make Ti and T2 in ljIII.4p and 
(|IIL5ll negligible, as the electrons move so fast that they can instantaneously adjust their motions with respect to the 
vibration and rotation movements of the slower and much heavier nuclei. 

The Born-Oppenheimer separation breaks down in several cases, chief amon g th em when the nuclear motion 
is strongly coupled to electronic motions, e.g., when the Jahn-Teller effects |43 are present. It also requires 
corrections for loosely held electrons such as those in Rydberg atoms. 

Research work on non-Born-Oppenheimer effects, the inclusion of the spin-orbit coupling, and on models without 
using the Born-Oppenheimer separation by treating the coupled dynamical motions of the electrons and nuclei 
simultaneously, may be found in Yarkony |48| and Ohrn [49|. for example. 



^^2Mk 2 \Rk-Rk'\ 

K,K' = 1 



B. Variational properties: the virial theorem and the Feynman— Hellman 

theorem 



The variational form of (IIII.lll is 



min *Lff * = E, 



(IIL7) 



where the trial wave functions ^' = 'i'{R,r) belong to a proper function space (which is actually the Sobolev space 
^i(]]j3(Afi+Af2)') ^j^g mathematical theory of partial differential equations (||21j. Chapter 2])). The (unique) solution 
^0 attaining the minimum of is called the ground state and the associated value Eq = (5'o|-ff|^o) is 

the corresponding ground state energy. Excited states ^'^ with successively higher energy levels may be obtained 
recursively through 



min(*|i?|*) = £;fc = 

where ^ is subject to the constraints 

(^-1^-) = 1, (^l^-o) = (^'l^'i) = • • • = Wk-i) - 0, for fc = 1, 2, 3, , 



(IIL8 



There are two useful theorems related to the above variational formulation: the virial theorem and the Feynman- 
Hellman theorem. We discuss them below. 

Let -0 be any trial wave function. The expectation value of the kinetic energy is 



fi2 



N2 



^1 2^^^ 



K 



^ 2m 



V. 



Now, consider the scaling of all spatial variables by 1 + A: 

7e=(-R,r)^(l + A)(-R,r) = (l 
Subject to the above transformation HIII.10|I . we have 

1 



\)n. 



(III.9) 



(III.IO) 



Ehi' 



(1 + A)^ 



By taking the variation of E^m with respect to A, we see that it is the same as taking the variation of E^m with 
respect to ^I*. Thus 



^5E. 



A=0 
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d 

dX 



(1 + A) 



2 ' 



A=0 



2(1 + A) • Ekin — —2Ekin — SEkin- (III.ll) 



A=0 



For the potential energy, the expectation is given by 

Epot = (*|^|*) = (V), (III.12) 



where V consists of the 3 summations of Coulomb potentials inside the bracket of ljIII.l|l (but V can be allowed to be 
any potential whose negative gradient is force). The spatial scaling (jllLlOp implies the change of displacement 6 = X, 
which is now regarded as infinitesimal. Therefore, 



= J ip*in)[n-vv{n)]tp{n)dn 



A=0 



R3" 



(7e • Vy) = Virial, (III. 13) 



where Virial is the quantum form of the classical virial. For ^E* ~ ^'o, the variational form l|III.7(l demands that 

S{H) = 6Ek,n + SEpot = 0. (III.14) 

Substituting HIII.11|) and HIII.13|) into (jIII.14p . we obtain 

~2EMn + Virial = 0. (III. 15) 

This is the general form of the virial theorem for an isolated system of an atom or a molecule. In the particular case 
that 

- "2 1 ^ 1 ^ ZkZk'c' 



fe,fe' = l K,K' = 1 

which is the power law 7?." (where TZ is the distance) with n = — 1, the classical virial property holds: 

Virial = nEpot = -E^ot (III. 16) 

as 

n ■ V7^" = {nTZ"-^)TZ ■ V7^ = ?l7^". 



^From piI.15|l and piI.16|l . we obtain the virial theorem for an (exact, not trial) wave function: 

2Ek.,n + Epot^0, (III.17) 

or 

Ekin = —-^Epot. (III. 18) 

This property is used as a check for accuracy of calculations and the properness of the choices of trial wave functions. 

The next, Feynman-Hellman theorem, shows how the energy of a system varies when the Hamiltonian changes. 

Assume that the Hamiltonian of an atom or molecule system depends on a parameter a, H — H{a). For example, 
a may represent the internuclear distance of the molecular ion Hj. The exact wave function 5* = ^'(q;) also depends 
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on a, so does the energy of the system E = E{a). Let's see how E{a) changes with respect to a, i.e., 

dE{a)/da^ — / -^i* {a)H{a)-^{a)dr 
da J 

= / — -^H{a)^{a)dr + / 'i'*(o;) ' ^{a)dr+ I ^*(a)H{a) ' dr 
J oa J oa J oa 

= E{a) I ^^^*(a)dr+ / m* {a)^^i^-^{a)dr + E{a) f ^*(a)^^^dr 
J oa J oa J oa 



E{a)-^{^{a)\-^ia)) + J **(c 
dH{a) 



da 
dH{a 



da 



'i'{a)dr 



*(a) 



da 



*(a) 



as (^'(a)l^'(Q:)) = 1 and is independent of a. The above is the Feynman-Hellman theorem. Its advantage is that 
oftentimes dH{a) / da is of a very simple form. For example, for the H^-like equation (jIV.13|) . upon taking a = R we 
have 



dH 
'dR 



and the average force on the nucleus B is 



d£ ^ ZgZb /Zbn\ 

dR i?3 \ r3 / 



That is the force and, hence, the potential energy curve requires calculation of (Zbr^/r^) only. This substantially 
simplifies the problem since the matrix element from the is no longer necessary. 



C. Fundamental properties of one and two-electron wave functions 

1. Riccati form, proximal and asymptotic conditions 

In this subsection, we introduce the Riccati form of the ground state wave functions as a unified way of understanding 
and deriving various cusp, asymptotic and correlation functions. This forms the basis by which compact wave functions 
for diatomic molecules can be derived. Consider the Schrodinger equation with a spherically symmetric potential in 
reduced units, 

-— V^VW + V(rU(r) = E^jM, (III.19) 

where /i = 1 is the central-force case, and fi — ^ is 

the equal-mass, relative coordinate case. We will be primary 
interested in 

Z 

I , (cf. pV.12p below) (III. 20) 

r 

however, the long range Coulomb potential is special in many ways and we can best understand the Coulomb-potential 
wave function by comparing and contrasting it to the short-range, Lennard- Jones potential 

y^j(-) = 'o(^-^]. (III.21) 



Since the ground state of (jIII.19|l is strictly positive and spherically symmetric, it can always be written as (unnor- 
malized) 



(III.22) 
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Substituting this into (|III.19p gives the Riccati equation for S{r) 



^V2^(r) - ^V5(r) • VS{r) + V{r) = Eq. (III.23) 

Since V5(r) = S"(r)f and V • f = 2/r, we have simply, 

-^5" + —S' - ^{S' f + Vir) = Eo. (III.24) 
2^ 2/i 

The advantage of this equation is that since the RHS is a constant Eq, all the singularities of V{r) must be cancelled by 
the derivatives of S{r). In particular, if we are only seeking an approximate ground state, then a reasonable criterion 
would be to require S(r) to cancel the most singular term in V{r). For both the Coulomb and the Lennard-Jones case, 
the potential is most singular as r — > 0. We will refer to this as the proximal limit. For both cases, the singularity of 
V{r) is only polynomial in r and therefore we can assume 

S{r) = ar", (III.25) 

giving explicitly 

— an(n + l)r"-2 _ J_^2^2^2n~2 ^ ^ ^ (III.26) 
2/i 2fi 

Note the structure of this equation, if the V{r) has a repulsive polynomial singularity, then it can always be cancelled 

by the term, leaving a less singular term r"~^ behind. For example, in the Lennard-Jones case, the most 
singular term is cancelled if we set 



— a Vr2"-2 + iL = 0, (III.27) 
2/i r^^ 

yielding, n — —5, and for fx — ^, the famous McMillan correlation function for quantum liquid helium. 



Sir) = IV^O 



r 



-5 



However, if Vir) has an attractive singularity, then it can only be cancelled by the r"~^ term, and if n < 0, would 
leave behind a more singular term instead. This means that an attractive potential V{r) cannot be as singular, or 
more singular, than — l/r'^, otherwise, the Schrodinger equation has no solutions. Fortunately, for the the Coulomb 
attraction (jIII.20p . the —Z/r singularity can be cancelled by setting 

—an(n + l)r"^^ - - = 0, (III.28) 
2/x r 

yielding, for /i = 1, n = 1, 

Sir) = Zr. 

For electron-electron repulsion with = ^ and Z = —1, we would have instead 

Sir) = -ir. 

In the Coulomb case, these proximal conditions are known as cusp conditions. A more thorough treatment of the cusp 
condition for the general case will be given in Subsection V.A and Appendices F and H. The proximal criterion of 
cancelling the leading singularity of Vir) can be applied generally to any Vir). The cusp conditions are just special 
cases for the Coulomb potential. Note that in ljIII.28p . the Coulomb singularity is actually cancelled only by the 
S '/ ipr) term of the Riccati equation, 

^S^-^^0, 
fir r 

since requiring S" to be a constant at the singular point forces 5" = 0. Thus the cusp condition can be stated most 
succinctly in term of the Riccati function Sir): wherever the nuclear charge is located, the radial derivative of Sir) 
at that point must be equal to the nuclear charge. 
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Next, we consider the asymptotic limit of r — > oo. In the case of short range potential, such that V{r) — > faster 
than 1/r, we can just completely ignore V{r) in (|111.24|l . Substituting in 



gives 



The constant terms determine 



S{r) ^ ar + /31n(r) 



-4-:if" + ^V + ^f- + ^H^o. (III.29) 



Setting the sum of 1/r terms zero gives, 



— a{l -(3)=0, (III.30) 
fir 



which fixes /3 = 1. The remaining terms will decay faster than 1/r and can be neglected in the large r limit. Thus 
the asymptotic wave function for any short-ranged potential (decays faster than 1 /r) must be of the form 

^(r) ie-"^ (III.31) 
r 

However, for the Coulomb potential —Z/r (we wish to reserve the possibility that Z can be distinct from Z and 
unrelated to Eq), we must retain it among the 1/r terms in (|111.30ll . 

— a(l - /3) - - = 0, (III.32) 
fir r 



resulting in 



/3=1-^. (III.33) 
a 



Thus the general asymptotic wave function for a Coulomb potential has a slower decay, 

^(0 ^ (III.34) 

For a single electron in a central Coulomb field, /i = l, Z = Z^a = Z,j3~Q^ and 

V'(r) ^ e-^^ (III.35) 

the decay is the slowest, very different from (IIIL31II . For more than one electron, or more than one nuclear charge, 
a ^ Z, (3 does not vanish and the correct asymptotic wave function is pil.34|) . We tend to forget this fact because 
we are too familiar with the single-electron wave function, which is the exception, rather than the norm. 

The proximal and asymptotic conditions are very stringent constraints: wave functions that can satisfy both are 
inevitably close to the exact wave functions. Satisfying the proximal condition alone is sufhcient to guarantee an 
excellent approximate ground state for all radial symmetric potentials such as the Lennard-Jones, the Yukawa {V = 
ie-"^Q^ > 0), and the Morse potential {V = e'^"'' - 2e-"'', a > 0). Needless to say, the proximal condition alone 
determines the exact ground state for the Coulomb and the harmonic oscillator potential. The significance of the 
proximal condition has always been recognized. The current interest |35l | in deriving compact wave functions for small 
atoms and molecules is based on a renewed appreciation of the importance of the correct asymptotics wave functions. 

2. The coalescence wave function 



Consider the case of two electrons orbiting a central Coulomb field 

1^2 '^^2 Z Z 1 

:5V? --V2 + — 

2 2 ri r2 ri2 



-l^l -1^1---- + —] V'(ri,r2) = E^Pirur2). (III.36) 
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Imagine that we assemble this atom one electron at a time. When we bring in the first electron, its energy is 
El — —Z'^/2, with wave function ip{ri) = exp(— Zri) localizing it near the origin. When electron 2 is still very far 
away, we can write the two-electron wave function as 

i/,(ri,r-2) = e-^''^"'^^''^^^ (111.37) 

Substituting this into (|III.36p yields the Riccati equation for S{r2), 

I5" - liS'f + ls'-- + — = Eo + \z\ (III.38) 

2 2 r2 r-i rx2 2 

The RHS defines the second electron's energy, E2 = Eq + ^Z'^ , whose magnitude is just the first removal or ionization 
energy. Since electron 2 is far away, and electron 1 is close to the origin, ri2 ~ r2- Thus electron 2 "sees" an effective 
Coulomb field —Z/r2 with Z = Z — 1. This is the case envisioned in ljIII.32|) with a — \J—2E2. The asymptotic wave 
function for the second electron is, therefore, 

^/'(ra) ^ r~''e-"''% (III.39) 

with /3 = 1 — (Z — l)/a. Since /3 is always less than one, the r^'' term is only a minor correction. Its effect can be 
accounted for by slightly altering a. The important point here is that the second electron need not have the same 
Coulomb wave function as the first electron. This coalescence scenario would suggest, after symmetrizing (|IIL37|I . the 
following two-electron wave function: 

^(n, rz) = e-^'^i-"'^^ -h e-"'-!-^'^^ (III.40) 

For the case where there are more than two electrons, one can imagine building up the atom or molecule sequentially 
one electron at a time. Each electron would then acquire a different Coulomb-potential wave function. This sequential, 
or coalescence scenario of approximating the ground state, in many cases, resulted in better wave functions than 
considering all the electrons simultaneously, which is the traditional Hartree-Fock point of view; see Subsection VI. C. 
In the case of He, the simple effective charge approximation 

with 

^cff = ^ - 4 = 1-6875, (III.41) 
16 

gives E^av — —Z^ff = —2.8476, while the "exact" value is -2.9037. The standard HF wave function of the form 

■0(^1,^2) = (f>iri)(f>{r2) 

improves [3^ the energy to Eyar = —2.8617. For comparison, the coalescent wave function Hill. 401) can achieve 
Eyar = —2.8674 at a = 1.286. Patil [s^, by restricting a to be consistent with the output variational energy via 
a = ^/—2Eyar — 4, obtained Eyai- = —2.8671 at a = 1.317. All these values of a are very close to the exact asymptotic 
value of a = \/— 2-^2 — •\/2(0.9037) = 1.344, lending credence to the coalescence construction. For arbitrary Z, by 
approximating Eg by —Z^g, we can estimate a by 



2Z2jj-Z2 i ^Z -5/8). (III.42) 

For Z = 2, this gives a = 1.30 (for small Z, we need to use the full expression rather than the approximation), an 
excellent estimate. This obviates the need for Patil's self-consistent procedure to determine a, and produces even 
slightly better results. The coalescent wave function l|III.40|) with this choice for a, defines a set of parameter-free 
two-electron wave functions for all Z. The resulting energy for Z = 2 — 10 is given in Table UTTl 
In 1930, Eckart has used wave functions of the form 

^(ri,r2) = e-"''i-'"'= +e-'"'i-'"'^ (III.43) 

to compute the energy of a two-electron Z-atom. His resulting energy functional is 



EEck{Z,a,b) = -Z{a + b) 



K{a, b) + ^(a^ + b^) + ab C(a, b) 



1 + C{a,b) [ ' ' ' 2 



(IH.44) 
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where 

a + (a + o)-' (a + by [a + o)° 

He obtained an energy minimum —2.8756613 for Re at a = 2.1832 and h = 1.1885. (We have used his energy 
expression to re-determine the energy minimum more accurately.) While the improvement in energy is a welcoming 
contribution, it seems difficult to interpret the resulting wave function physically. How can an electron "sees" a 
nucleus with charge greater than 2? 

To gain further insight into Eckart's result, and coalescence wave function in general, we note that pil.43|l can be 
rewritten as 

X{ri,r2) = e-^(''i+''^) cosh [B(ri - ra)] , (ni.46) 

with A = (a + 6)/2 and B = (a- b)l2. This form has the HF part e-'^(''i+''2). If we substitute Eckart's values, we 
see that A — (2.1832+ 1.1885)/2 — 1.6859, which is nearly identical to the effective charge value (|HI.41|) . This is not 
accidental. If we use the approximation (|III.42p for a, then the coalescence wave function (|III.40ll would automatically 
predict 

A=\{Z + a) = Z-^\ (III.47) 

Thus within the class of Eckart wave function 1111.431) , the coalescence scenario correctly predicts the path of optimal 
energy as being along A = Z^ff. Moreover, this improvement in energy, which has laid dormant in Eckart's result 
for three quarters of a century, can now be understood as due to the radial correlation cosh term in ljIII.46|) , built-in 
automatically by the coalescence construction. This term is the smallest (=1) when ri = r2, but is large when the 
separation ri — r2 is large, i.e., it encourages the two electrons to be separated in the radial direction. 

This suggests that we should reexamine Eckart's energy functional in terms of parameters A and B. Expanding 
(|III.44p to fourth order in B yields, 

EEckiz, A, B) = -zls + {A- z,sf -\y + - (ni.48) 

where y = B"^ /A. li B = A = Z^s this yields the effective charge energy —Z'^q. Regarding the effect of B as 
perturbing on this fixed choice of A, the 1/A « 1/Z term can first be ignored. Minimizing y simply yields 

Ee..{z,a,b) = -z..-^^-^^ im.49) 

at 

where we have restored the 1/A term. This remarkably simple result is the content of Eckart's energy functional. The 
energy is lower from the effective charge value by a nearly constant amount 3/128. In column three of Table ITTll we 
compare this approximate energy l|III.49p . with the absolute minimum of the Eckart's energy functional on the fourth 
column. The agreement is uniformly excellent. By comparison, we see that the coalescence construction, wz</ioMi 
invoking any minimization process, also gives a very good account of the energy minimum. 



3. Electron correlation functions 



Coalescence wave functions are better than HF wave function because they have built-in radial correlations. To 
further improve our description of He, as first realized by Hylleraas '30], one can introduce electron-electron correlation 
directly by forcing the two-electron wave function to depend on ri2 explicitly. (A more detailed discussion of correlation 
fmictions will be given in Subsection V.B below.) Again, our analysis is simplified by use of the Riccati function. Let 



V'(ri,r2) =e-^('-i'^^''^^=\ 



(IH.51) 
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TABLE III: Ground state energy of the He-like ions as calculated from various wave functions. 
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but now consider 

5(ri,r2,ri2) =Zri+Zr2+.9(ri2). (III.52) 

We have 

Vi5 = Zvi+ g'r 12,^28 = Zv2+ g'v2i, 
^2^ _ 2Z , 2.9' , 2Z , 2g' 



^IS = — + -^+g",WlS^— + ^+g\ 
n ri2 r2 ri2 



and H111.36() in terms of 5* reads 



9" + ^^-^ - {gT - Zg'ih - V2) ■ fi2 -Z^ = Eo. (III.53) 

ri2 

In order to eliminate the 1 /ri2 singularity, we must have 

lim g'{r 12) = -i. 
1-12^0 2 

Thus, one can consider a series expansion for g{r 12) starting out as 

9(ri2) = -^ri2 + icr^2 + ••• 
Keeping only up to the quadratic term, in the limit of ri2 — > 0, (|III.53|I reads 

3C - ^ - Z^ + Oiri2) = Eo. 

If <?(''i2) were exact, the LHS above would be the constant ground state energy for all values of ri2. Inverting the 
argument, we can exploit this fact to determine C at ri2 = 0, provided that we can estimate the ground state energy 
Eq. The simplest estimate for a two electron atom would be Eq = — Z^, implying that 

C=^. (III.54) 

However, since the effective charge approximation for the energy is much better, we should take instead, Eq = —Z^^, 
thus fixing 

C^^ + i(Z=-Zi,).l + |(z-|)^ (....55) 

The determination of the quadratic term of g{r 12) was advanced only recently by Kleinekathofer et al. [sJl- (If we also 
improve the one-electron wave function from cxp(— Zr) — > exp(— Zcffr), then C must go back to the value (|III.54p . 
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The coefficient C therefore depends on the quahty of the single electron wave function.) The wave function ljIII.51|) 
can now be written as 



V'(ri,r2) =e-^'^^e-^^V(^i2) 



where 



/(ria) = exp 



-ri2(l - Cria) 



(III. 56) 



Since the above argument is only valid for small ri2, the large ri2 behavior of f(ri2) is not determined. It seems 
reasonable, however, barring any long range Coulomb effect, to assume 



lim /(ri2) 



constant. 

The form l|III.56l) can have behavior l|III.57|) if we just rewrite it as 

ri2 



fpj{r 12) = exp 



2(1 + Cri2) 



This Pade-Jastrow form has been used extensively in Monte Carlo calculations of atomic systems 
to achieve (|III.57ll , Patil's group [s^ |^ have suggested the form 



(III.57) 

(III.58) 
Alternatively, 



fp{ri2) = — (l + 2A-e-^'-i^); cf. 06.(3)) andFig.[I3 



which has the small ri2 expansion 



fp{ri2) = 1 + -ri2 



' 12 



By comparing this with similar expansion of (jIII.56p . we can identify 



X = 2C- 



, Z 

12 V 32 



Le Sech's group [13 have employed the form 



with 



(III.59) 



(III.60) 



(III.61) 



(III.62) 



a — -A. 
2 

This function is not monotonic; it reaches a maximum at ri2 = 1/a before level off back to unity; cf. Fig. 1111 in 
Subsection 0B. However, this point may not be practically relevant, since most electron separations do not reach 
beyond the maximum. For the sake of comparison, we may use the maximum of Le Sech's correlation function as 
its asymptotic limit. This is the natural thing to do because all three functions can now be characterized by their 
asymptotic value as ri2 00 : 



f p. j{r 12) exp 



1 

2C 



1 



A + 1/2' 



/p(^i2) 
fhsir 12) 



1+2A' 
1 

1 + ^- 



Their approaches toward unity are approximately A"^, ^A^^, and ^A~^, respectively. Note also that as A increases 
with Z , the asymptotic value of f{ri2) decreases, this is the correct trend long observed in Monte Carlo calculations 
on atomic systems [S^l- For Z = 2,we take C = 1/2, A = 1/2, a = 1/4 and compare all three correlation functions in 
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Pade-Jastrow 




FIG. 7: Comparison of three electron-electron correlation functions. 



Fig. Also plotted is the simple linear and quadratic forms. The simplest linear correlation hmction, (1 + ^''12), is 
very distinct from the other three. 

We now have all the ingredients needed to construct an optimal wave function for He. First, with the introduction 
of explicit electron correlation function, there is no need for the radial correlation introduced by the coalescence wave 
function, i.e., we are free to abandon the radial correlation function cosh(_B(ri — r2)). Second, the asymptotic form 
of the wave function exp(— ar) should be maintained. However, this asymptotic wave function, when extended back 
to small r, violates the cusp condition. These concerns can be simultaneously alleviated by replacing the second 
electron's wave function via 



-Zri 



cosh(/3r2). 



At small r2, the cosh function is second order in r-i and therefore will not affect the cusp condition. 
cosh(/3r2) exp(/3r2), and the choice 



At large r-i-, 



(ni.63) 



will give back the correct asymptotic wave function. Upon symmetrization, we finally arrived at the following compact 
wave function for He: 



-Zr 



^'^^^ [cosh(/3ri) + cosh(/3r2)] /(ri2). 



(HI.64) 



This wave function, first derived by Le Sech [5j|, satisfies all the cusp and asymptotic conditions. We fixed a, /3 and 
C by ljHI.42|) . IIHI.63|) and (jIII.55|l . respectively, and there are no free parameters. The only arbitrariness is the form 
the correlation function /(ri2). Since fp{r 12) is bracketed by fpj{r 12) and /Ls(ri2), we only need to consider the 
latter two cases. For Z — 2 — 10, the resulting ground state energy for the two electron atoms are given in Table ITTll 
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4- The one-electron homonuclear wave function 



The one-electron, homonuclear two-center Schrodinger equation 

1 V2 - - - -) V^r) = i?^(r), (III.65) 

where Va + R./2|, r?, = |r — R/2|, can be solved exactly, as will be shown in Subsection IIVI B below. However, its 
ground state wave function can also be accurately prescribed by proximal and asymptotic conditions. For Z = 1, this 
is hydrogen molecular ion problem. As first pointed out by Guillemin and Zener (GZ) [54 |. when i? = 0, the exact 
wave function is 

^(r) = e"2^'' (III.66) 

and when R ^ oo, the exact wave function is 

^(r) = e'^'-" 4- e-^'■^ (III.67) 

They, therefore, propose a wave function that can interpolate between the two, 

i'Gz{ra,ra) = e-^i'-^-^^'-" + e-^^'■»-^l'^^ (III.68) 

At i? = 0, Tf, = rf, = r, one can take Z\ — Zi — Z . At i? — > oo, one can choose Z\ = Z and Zi = 0. At intermediate 
values of i?, Z\ and Z2 can be determined variationally. This GZ wave function gives an excellent description [s^ 
of the ground state of HJ. As explained by Patil, Tang and Toennies another reason why this is a good wave 
function is that that ljIII.68(l can satisfy both the cusp and the asymptotic condition. We can simplify Patil tt al's 
discussion by rewriting the GZ wave function, again, in the form 

i^Gz{ra,n) = e-^('-"+'^'')/2 cosh [B{ra ~ n)/2] , (III.69) 

when R ^ 0, A ~ 2Z, and when R —f 00, A = B = Z . The imposition of the cusp condition can be done most easily 
in terms of the Riccati fmiction. We, therefore, write 

^Gz(ra,rb) = e-^('--'-'), (III.70) 

with 

S{ra,n) = Aira + rfc)/2 - ln( cosh[B(r, - n)/2] ). 
The cusp condition at = is then easily computed, 

dS 



= Z, A + B tanh(Bi?/2) = 2Z. (III.71) 

ra=0,rb=R. 



One can verify that this is also the cusp condition at ri, = 0. From ljIII.7ip . one sees easily that a.t R — 0, A — 2Z, 
and when R ^ 00, A + B = 2Z. At finite R, the asymptotic limit r 00 means that = r;, = r and the GZ wave 
function approaches 

^Gz{ra,n) -> e"^*^. 
On the other hand, the exact wave function must be of the form l|III.34ll 



^{ra, n) ^ e-V2^E^--/3 !"('■), (III.72) 

with (3 — 1 — 2ZI ^2\Eq\. Since 2Z > y/2\Ea\ > Z, we can estimate that at intermediate values of R, y/2\Ea\ ~ ^Z, 
suggesting a negligible /3 ~ — i. Thus it is suffice to take 



A « y^2\Eo\ (III.73) 

Guillemin and Zener have allowed both A and B to be variational parameters. Patil et aPs estimate of A is 
essentially that of l|III.73l) but with slight improvement to incorporate the variation due to /3 0. We adhere to the 
cusp condition 1)111. 71|l but allow A to vary. In practice, it is easier to just let B vary and fix A via the cusp condition 
(ITTT.7H . In all cases, this wave function can provide an excellent description of the hydrogen molecular ion, with 
energy derivation only on the order of 10"'^ Hartree over the range of i? = — 5. 
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5. The two-electron homonuclear wave function 



The two-electron homonuclear Schrodinger equation is, 



^V? -iv^- — - — - + v.(ri,r2) = i?V(ri,r2). (III.74) 

2 2 ria rit ri2 / 

Let's denote the one-electron two-center GZ wave function as 

0(r) = e^'^'^ cosh(B,5) 

where we have defined 

= 7^ 5 = . 

2 2 

(The variables ti and 5 here will correspond, respectively, to A and /i of the prolate spheroidal coordinates in Subsection 
IIVI B.) To describe the two-electron wave function, if one were to follow the usual approach, one would begin by defining 
the Hartree-Fock like wave function 

7/'(ri, rs) = 0(ri)0(r2) = e-^(''i+'"=) cosh(Mi) cosh(S(52). (III.75) 

However, this molecular orbital approach is well known not to give the correct dissociation limit of H2. In the limit 
of i? ^ CX3, we know that 

-Zra I „-Zrt 



g-Zr„ _^ g-Zr. (III.76) 



and therefore 



7/.(ri,r2) ^ (e'^'-i" +e-^'-")(e-^'^^'' -f e-^''^") 

^ (e'^'^i-e-^''^'' + e-^'-i''e-^''=^°) -I- (e-^'-i-e"^''^" + e'^^'^^e-^''^^) (III.77) 

Only the first parenthesis, the Heilter-London wave function, gives the correct energy of two well separated atoms with 
energy E — 2{—^Z^). The remaining parenthesis describes, in the case of H2, the ionic configuration of H^H+, which 
has higher energy than two separated neutral hydrogen atoms. Thus the molecular orbital approach ljIII.75|) will alway 
overshoot the correct dissociation limit. This is a fundamental shortcoming of the molecular orbital approach and 
cannot be cured by merely improving the one-electron wave function, i.e., by use of the exact one-electron, two-center 
wave function. Even the coalescent construction cannot overcome this fundamental problem. In the large R limit, 
the inner electron's wave function must be ljIII.76p . and hence no matter how one constructs the outer electron's 
asymptotic wave function, one can never reproduces the Heilter-London wave function. In both the molecular orbital 
and the coalescent approach, one must resort to configuration interaction to achieve the correct dissociation limit. 
Even if one were to use the exact one-electron wave function in doing configuration intereaction, as it was done by 
Siebbeles and Le Sech jl^, the energy still overshoots the correct dissociation limit if the correlation (1 -t- ^ri2) is 
used. This is because we must have /(ri2) — + const in order to reproduce the Heilter-London limit. 

However, one can learn from Guillemin and Zener's approach, and insist on a wave function that is correct in both 
the R = and R ^ 00 limit. This seemed a very stringent requirement, but surprisingly, it is possible. The wave 
function is 

^(n, rs) = 6-^(^1+'^^) cosh [B{Si ~ 62)] /(ri2). (IH.78) 

For R — 0, ai — ri, a2 — r2, and the above function reduces to 

V'(ri,r2) =e-^('^i+'-^)/(ri2). 

which is not a bad description of He. In the limit of i? ^ 00, if we take A = B = Z , we have 

^(ri,r2) = (e-^'^i-e-^'-^^ +e-^''"e-^'^^")/(ri2), 

= (e-^'-i-e"^''^^ -Ke-^''"e-^''^°), (IH.79) 

since f{oo) — > 1. Thus wave function (|III.78p is the simplest homonuclear two electron wave function that can describe 
both limits adequately. The wave function ljIII.78|) for H2 without /(ri2) has been derived some time ago by Inui [s^ 
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and Nordsieck jSq . However, they were only interested in improving the wave function and energy at the equilibrium 
separation and were not concerned with whether the wave function can yield the correct dissociation limit. 

To estimate the form of the correlation function /(ri2), we repeat our analysis as in the Helium case. The two- 
electron wave function can again be written in the Riccati form (|in.51ll , but now with 

5(ri, r2, ri2) = Zri„ + Zr^b + .g(ri2), (HI.SO) 

where we have assumed the unsymmetrized form of the Heilter-London wave function. The resulting equation for g 
is also similar to (|in.53ll . 



/ + {g'f - Zg'iha - r2b) ■h2-Z^ + 0(-] =Eo. (IH.81) 

ri2 \RJ 

In the Helium case, the dot product term vanishes in the hmit of ri2 — > 0, here it does not. In the case where the two 
electrons meet along the molecular axis, f = z, f2fc = — z, fi2 — — z, the resulting equation 

V + 1 , ,s2 , ' v2 , 



+ (gf + 2Zg' ~Z^ + 0{-] =E„ (ni.82) 

ri2 \R; 



can be solved by setting Eq — —Z\ and expanding 



In the hmit of ri2 0, Hill. 82(1 reads 



.9(^2) = - \t\-2. + 



3C-i-Z + O(ri2) = 0, 



giving, 



This agrees with Patil et al's result 16] of 



+ (ni.83) 



A = 2C-i = i(2Z-l), 

but without the need of consulting hypergeometric functions. For the H2 case, C = 5/12 = 0.42. In our calculation 
with wave function ljIII.78|) . with fpj{r 12) given by ljIII.58|) . the energy minimum at intermediate values of R is at 
C = 0.40, in excellent agreement with the predicted value. Since C « 0.50 for Helium, C's variation with R is very 
mild. 

The resulting energy for the wave function l|III.78|) . is given in Fig.|Sl(solid line). We vary the parameter B, while 
the other parameters A and C are fixed by Eqs. Ijlll.7ip and l|III.83|) . The parameter B is 0.8 for R < 2, and 
moves graduately toward one at larger values of R. The energy at equilibrium is as good as Siebbeles and Le Sech's 
calculation with unsealed HJ wave functions and correlation function (l-l-^ri2) (triangles). Without configuration 
interaction, Siebbeles and Le Sech's energy overshot the dissociatin limit as shown. The wave function ljIII.78|) can 
be further improved by adding a coalescense component a la Patil, Tang and Toennies 0|. This will be detailed in 
the next subsection. 



6. Construction of trial wave functions by Patil and coworkers 



The previous discussion demonstrates that relatively simple wave functions which incorporate the cusp conditions 
and the large distance asymptotics, having no or only a few variational parameters, can be constructed to yield fairly 
accurate results. Here we mention other similar trial wave functions studied in the literature. The importance of the 
local properties in the calculation of the chemical bond has been emphasized by Patil and coworkers |16l 1371 13M I3M liol 
IH^ . Their analysis is in the spirit of our previous discussion, however, yields more complicated trial wave functions. 
Here we briefly mention the main aspects of the construction scheme and, to be specific, consider the ground state of 
H2 molecule described by the Schrodinger equation (IIII.74|I . 
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FIG. 8: Ground state potential curve E{R) of H2 for different trial wave functions. Triangles correspond to "Le Sech" molecular 
orbital calculation with exact one-electron orbitals. Our wave function 1111.7811 (solid line) gives the correct dissociation 
limit. Patil et al.'s results are shown as small dot (function 1111. SSL 1111. 91L 1111.9211 ') and dash (function 1111.1021 1') lines. Large 
dots are "exact" values of Ref. 0- 



Let us assume that r2 3> ^i, i?, 1. Then ri2 ~ r2a ~ ^26 ~ ''2 and the HamUtonian reads 



H 



1 



z 



Z , 1^, (2^-1) ^ 

ria rib R 2 ^ r2 ' 



(in.84) 



The first four terms in (|in.84|) yield the problem, while the last two terms correspond to motion of a particle 
in a Coulomb potential with an effective charge 2Z — 1. This Hamiltonian allows us to separate variables and write 
\l/(ri,r2) = (ri)i^(r2), where the function (/j(r2) satisfies the equation 



1 o (2Z- 1) , , , 



(in.85) 



e = £'jj+ — i? > is the ionization energy of the H2 molecule, and E is the ground state energy of H2. As a result, 
the asymptotic behavior of \E'(ri,r2) at r2 3> ri, R, 1, is 



^'(ri,r2) 



(2Z-1)/V2i-1 



exp(-V2£r2)*H+(i"i), 



(in.86) 



which is similar to the coalescence wave function (IIIL39I) for He. Now assume that r2 ^ ri 3> -R, 1, then 5'jj+(ri) is 
given by Eq. HIII.72P and, therefore. 



,T,/ ^ 2Z/V2£r-1 (2Z-1)/V2s-1 I /K /TP N 

*(ri,r2)wri r\ " exp(-V2£iri - V2er2), 



(in.87) 
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where Si = Z'^/R — > is the separation energy of electron in Hj. The power-law factor slowly varies as 
compared to the exponential decaying contribution. Hence, one can assume the power-law factor to be a constant or 
approximate the combination r" exp(— 5r) as 

r" exp(— for) = exp(— 6r + alnr) w exp(— for + alnr'o + (r — ro)a/ro), 

where rg can be determined as a variational parameter or chosen to be ro = i? -I- 1/6 [l^ . 

To incorporate the cusp conditions and the large distance asymptotic the trial wave function is separated into two 
parts 

*(ri,r2) = $(ri,r2)/p(ri2), (III.88) 

where fp{r 12) is the Patil et al. electron-electron correlation function given by Eq. (|III.59|I . Roothaan and Weiss js^ 
have made a very accurate numerical investigation of the desired correlation function for the ground state of the He 
atom. In the vicinity of ri2 = 0, the correlation function is linear and satisfies the cusp condition. It monotonically 
increases and approaches a constant as ri2 becomes very large. Clearly the function fp{r 12) satisfies these conditions 
(see Fig. I12|l . In the united atom limit {R = 0) it was found that the energies computed with the variationally 
determined A are essentially the same as given by the analytical expression, 

A = 4 (IH.89) 

derived from a theory in which l/ri2 is treated as a perturbation 37]. In a molecular system, as R increases, one 
should expect A to decrease monotonically and become vanishingly small for R 00. The small and large R behavior 
is satisfied provided [s^ 

A^ ^^f-y^ , (III.90) 

l + 10Z3i?2/(15Z-6) ^ ' 

The electron- nucleus cusp conditions do not uniquely define the space wave function $(ri,r2). If one wishes to 
maintain the electronic configuration idea with an independent particle picture, one can adopt the following form of 
^'(ri,r2): 

$(ri,r2) =(/.(ri)0(r2), (III.91) 

with <t>{rj), i — 1, 2, being the Guillemin-Zener |54| trial wave function for H^ 

0(rj) = exp(-2:irja - Z2rjb) + exp{~Z2rja - Zirjh), (III.92) 

where 21 > 0, Z2 > are variational parameters. Alternatively, zi and Z2 can be determined from the cusp conditions 
at rja — and r^b = for j = 1, 2. The wave function ljIII.9ip . (|III.92p is identical to ljIII.75|) which is known not to 
give the correct dissociation limit of H2. However, for the pedagogical reason we briefiy discuss it here. 
As ria approaches zero 

$(ri,r2) ^ 0(r2) [(1 - ziri,,) exp(-Z2i?) + (1 - Z2ria) exp(-2ii?)] . (III.93) 

Imposing the cusp condition ^(ri, r2) — > G'(r2)(l — Zria) we obtain an equation for zi and Z2'- 

zi = Z+{Z- Z2) exp[~(zi - Z2)R]. (IH.94) 

Thus, if zi and Z2 are related as in Eq. Ijlll.94p . then the electron- nucleus cusp conditions are automatically satisfied. 
The second equation for zi and Z2 can be determined by the asymptotic condition. For r2 3> 3> i?, 1 Eqs. (|III.9ip . 
(Iin.92^ yield 

*(ri, r2) w exp[-(zi -f Z2)ri - (zi -I- Z2)r2]. (III.95) 

^From the other hand, according to Eq. (|111.87p . the wave function must have the following exponential behavior 
\I'(ri,r2) ~ exp(— Y^2e7ri — \/2er2)- The two parameters zi and Z2 do not allow to match the asymptotic exactly. 
However, one can approximately choose (34| 

^2 

zi + Z2=ei+e= — -E. (in.96) 
H 
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Eqs. Ijlll.94|) . Ijlll.96|) determine zi and Z2 self-consistently together with the ground state energy E. 

Kleinckathofcr ct ah 34] used the trial function H111.88() . (|111.91|l . Hill. 92(1 with A, zi and Z2 determined from Eqs. 
(iTTmi^ . fTM . mm . The wave function has no free parameters and yields 4.661 eV for the binding energy of the 
H2 molecule which is very close to the exact value of 4.745 eV. However, E{R) becomes less accurate at large R and 
fails to describe the dissociation limit. The corresponding E{R) is shown as a small dot line (Patil et al. 1) in Fig. |S1 

Both the cusp conditions and the large distance asymptotic can be satisfied exac tly p rovided more sophisticated trial 
functions are introduced. For example, for small and intermediate R, Patil et al. suggested to use a combination 
of "inner" and "outer" molecular orbitals which are build from the Guillemin-Zener one-electron wave functions: 



^'„(ri,r2) = [<?!)i„(ri)(/)out(r2) + (r2)<?!>out(ri)]/p(7-i2), (III.97) 

where the "inner" orbital is 

(^i„(rj) = exp(-zirja - Z2rjb) + exp(-Z2?'ia - zirjb)- (III. 98) 

Analogously, an "outer" orbital is defined as 

(/>out(rj) = exp(-Z3rja - z^rjb) + cxp(-Z4rja - zsrjb). (III.99) 

All the parameters zi, Z2, Z3 and Z4 are determined by the cusp and asymptotic conditions. 

At large R, the atomic orbital wave function provides a better description of the two electron system. The appro- 
priate wave function is 

^-alri, r2) - [$(ri, r2) + $(r2, ri)]/p(ri2), (Ill.lOO) 

where 



$(ri,r2) = exp[-Z(ria + ru + r2a + r2b)] {cosh(z5ri6) cosh(z6r2Q) + cosh(z6ri(,) cosh(z5r2Q)} . (III. 101) 

Eq. ljIII.100|) satisfies all the electron-nucleus cusp conditions. At the same time, it has two free parameters Z5 and 
Ze which can be used to satisfy the two asymptotic conditions. 

For a description in the entire range of internuclear distances, one can use a linear combination of the two wave 
hmctions just discussed 

* = *™-M)*a, (III.102) 

where D is a variational parameter. For II2 the molecular orbital '^m dominates in the region R < 1.7, while the 
atomic orbital 'i'a dominates at i? > 1.7. With this complicated one parameter wave function Patil et al. |1_6] obtained 
4.716 eV for the binding energy of H2 molecule and a very accurate potential curve in the entire range of R. The 
corresponding E{R) is shown as a dash line (Patil et al. 2) in Fig. |S1 Similar wave functions which take full advantage 
of the asymptotic and proximal boundary conditions are useful in variational calculations of larger systems |39l | . 



IV. ANALYTICAL WAVE MECHANICAL SOLUTIONS FOR ONE ELECTRON MOLECULES 



From now on throughout the Sections IIVI IVI and IVll unless otherwise noted, we assume the Born-Oppenheimer 
separation, where there are N nuclei, containing Zk protons located at Rk, respectively, for k — 1,2, . . . , N, and 
electrons. Each electron's coordinates are denoted as rj, j = 1,2, .. . ,Ne, where rj = {xj,yj, zj). The steady-state 
equation, in atomic units, can be written as 

j=l l<j<k<N^ ■' i=l fe=l ' ■> ' l<j<fc<7V ' ■> ' 

where 

q2 q2 q2 

^/> = V(ri,r2,...,rwJ, = «^ + ^ + rok = \r,-rk\, j,k = 1,2, . . . ,N^. 



We wish to solve the eigenvalue problem (jIV.ip . 

Closed-form solutions to (jIV.ip are hard to come by in general. What is known today is the following: 
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(i) JVe = l,JV=l: 

This is the case of hydrogen atom H (or H-hke ions with a single nucleus and electron) whose solutions are 
known explicitly in closed form, to be briefly reviewed in Subsection IIV Al 

(ii) Ne = l,N = 2: 

This is the (or H^-like) two-centered molecular ion, whose solutions are separable and expressible as an 
infinite series of products of special functions in prolate spheroidal coordinates, where coefficients of the series 
are not explicitly given. Such are the two renowned classic solutions due to Hylleraas and Jaffe, to be discussed 
m Subsection IrVBl 

(iii) Ne = l,N>3: 

This one-electron multi-centered molecular ion has an analytic solution due to Shibuya and Wulfman |59| | in 
terms of integral equations on the unit hypersphere of the 4-dimcnsional momentum space. This will be reviewed 
m Subsection ITVOI 

Except for Case (i) above, one must resort to numerical methods in order to derive quantitative and qualitative 
information, for all cases where ^ 1, N ^ I. Our particular interest in this paper is the diatomic case, with 
= N = 2, using the orbitals in Cases (i) and (ii) above as the building blocks. 

A. The hydrogen atom 

When Ne = N = 1, with Zi = 1 and Ri = 0, equation (|IV.1|I becomes 

-iv^ - -] V'(r) = Ei;(r), r e M^ (IV.2) 
2 r J 

the Born-Oppenheimer separation of the hydrogen atom. 

We write (IIV.2|) in spherical coordinates in view of the symmetry involved: 



where 



" ^ ^OYe G^^'l) ' ^^^^ Legendrian); (IV.4) 

r = the radial variable, < r < oo; 
6 — the colatitude, < < tt; 
6 — the azimuth, < (t> < 27r. 



Equation IIIV.2|) has separable solutions 

^{r) ^^{r,e,(j)) = R{r)Y{e,(j3). (IV.5) 
The angular variables are quantized first as we know that angular functions are the spherical harmonics 

¥{9, (j)) = Yirnie, 0) = e,™(0)$™(0), e = 0,l, 2,. ..,m = e,£-l,. (IV.6) 
on the unit sphere ^2 = {r G R'^ | = 1}, satisfying 

A2y,„(^^, 0) = + i)Yim{e, </)), (iv.7) 

where in l|IV.6|l . 

$„(0) = (2^)-^e™^ (IV.8) 

f (2£+l)(£-|m|)! )^ j^l^ 
Qem[d) ^ < , I |.| > Pi '(cos 6*), (IV.9) 



2(^+|m|)! 

(pj™' is the associated Legendre function). 
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Using ljIV.5|l - (|IV.7|l in ljIV.2|l . we obtain the equation for the radial function 



\(^rR)" -^-^:^^R+(^~2E]R = 0. (IV.IO) 



Solutions to the eigenvalue problem (|IV.10|I that are square integrable over < r < oo are known to be 

Rir) = R,Ar) = (-2) {^^[(^^^j (2r)^L^t+/(20e-^ (IV.U) 

En — n — 1,2, ... , independent of i, 

2 

where L^^+i ^''^^ ^"^^ associated Laguerre functions such that for to, ?i = 0, 1, 2, . . . 

xLZ+n" + (m + 1 - x)LZ+r^' + i^n + n)LZ+,, = 0, 

pX —(m+n) jm+n 

™+"^ ^" (TO + n)! 

(when TO = 0, L^ix) is simply denoted as Ln{x)). 

In the subsequent sections, we will utilize mainly the ground, or Is state of the hydrogen atom, where n — 1,1 — Q, 
i.e., 

$(r) = ^e"''; cf. Fig. HOI (IV.12) 
v2 

B. /fj^-like molecular ion in prolate spheroidal coordinates 

We now consider the eigenvalue problem for two-centered H^-like molecular ion with one electron and two fixed 
nuclei with effective charges Za and Z},. Given R the internuclear separation distance, we want to find E and ^E* such 
that 

^(^ + ^- ^] VP = E^. (IV.13) 

2 \ra n R J 

In Appendix A we show how to separate the variables through the use of the ellipsoidal (or, prolate spheroidal) 
coordinates (see Fig. E)) 

^=1 v/(A2-1)(1-m2) COS0, y=|V(A2-I)(l-^2) z = |am. (IV.14) 

In such coordinates the wave function can be written as 

* = A(A)M(/x)e*™'^. (IV.15) 

Separation of variables yields 



^{(A^-l)^} + {A + 2i?,A-,^A^-3^|A = 0, i?r.^^^^, A>1; (IV.16) 



' ^l-M^)^U(-^-2i?,M + pV-T^|M = 0, i?. . ^^i^V^, iMl < 1. (IV.17) 



dfi y "/^ J L 1 — /i"^ 

Note that A and p are unknown and must be solved from ljIV.161) and (|IV.17|) as eigenvalues of the coupled system. 
Once A and p are solved, then E can be obtained from (|A.8|I . 

In the next two subsections, we address the issues of solving (jIV.16|l and (|IV.17|I . respectively. 
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FIG. 9: (a) Elliptical coordinates (A, /i). (b) Prolate spheroidal coordinates (A, /i, (p) with A = {ra + r%)/ R and /i = [fa — rt,)/R. 
The range of coordinates is 1 < A < oo, —1 < < 1 and < < 2-k. 

1. Solution of the K- equation 11V.16I I 

To solve (jIV.16|l . it is important to understand the asymptotics of the solution. Rewrite l|IV.16|l as 

(A^ - 1)A"(A) + 2AA'(A) +(a + 2Ri\ - p^X^ - j^^) A(A) = 0. (IV.18) 



First, consider the case A ^ 1; we have 

^^ + 2^^^ - 1a(A) (IV.19) 



2A 

= A"(A) + ^^A'(A) + 



A2 - 1 ^ A2 - 1 (A2 - 1)2 

« A"(A)-/A(A), for A > 1. 

This gives 

A(A) « aie-P^ + oze^^, for A > 1, where p > 0. (IV.20) 
But the term 026^^ has exponential growth for large A, which is physically inadmissible and must be discarded. Thus 

A(A) « aie-P^, for A > 1. (IV.21) 
A finer estimate than ljIV.21|l can be stated as follows 

00 

A(A) = aoe-P^A^ J^TJ » 1) (^^.22) 

j=o 

where (3 = Ri/p— 1, cq = 1, ci = (p2 — f]'^ — 2/3)/(2i?i — 2p(3 — p), and oq is an arbitrary constant. Proof is given in 
Appendix ^ 

Next, we consider the case A > 1 but A ~ 1. In such a limit we have (see Appendix O 



00 



A(A)«(A-I)^^c,(A-l)'=. (IV.23) 



k=0 
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Our results in ljIV.2ip and ljIV.23|l suggest that the form 

A(A) = e~P^{X - l)^Xl^f{X), for some function /(A) (IV.24) 
would contain the right asymptotics for both A ^ 1 and Awl. Here, obviously, /(A) must satisfy 

/(I) ^ 0, lim |/(A)| < C, for some constant C > 0. (IV.25) 



A — >oo 



Actually, in the literature ([H0)I^E3)) improved or variant forms of the substitution of IIIV.24p are found to be 
most useful: 

(i) (Jaffe's solution |6l|') 

A(A) = e-P\\^ - 1)^ (A + 1)'^ £ 5n (xTl)" ' ('^ ^ y - - l) ■ (IV.26) 

n— ^ \ y / 

This leads to a 3-term recurrence relation 

anQn-l - Pngn + 1ngn+l = 0; 71 = 0,1,2,...; 9-1=0, (IV.27) 

where 

q;„ = (n — 1 — (T)(n — 1 — (7 — to), 

^„ = 2n2 + (4p- 2cr)n- _2pcr- (to + l)(TO + cr), (IV.28) 

7n = + + m + 1), 



and, consequently, the continued fraction 



70 p. 



^ 72^3 
P2 - ^5 



for A and p. 

(ii) (Hylleraas' solution [30]) 



oo 

A(A) = e-P^'-'Hx' - 1)^ 5] .^^^L^^^i^), X EE 2p{\ - 1), (IV.30) 

•■^ — ' TO + n ! 

n— ^ ' 

where i™+„ is the associated Laguerre polynomial and c„ satisfy the 3-term recurrence relation 

a„c„_i - /3„c„ + 7„c„+i = 0, n = 0,1,2,...; c_i = 0, (IV.31) 

where 

O-n = — Tn) (n — TO — 1 — cr) , 1 

/3„ = 2(n - to)2 + 2(n - to)(2p -a)-[A-p^ + 2pa + (m + l)(m + a)], > (IV.32) 
-Jn = {n + l){n - 2m - a), J 

and the same form of continued fractions ljIV.29|) . 

2. Solution of the M-equation 11 V. 1711 

Equation (|IV.17|) has close resemblance in form with (|IV.16|) and, thus, it can almost be expected that the way to 
solve (jlV.iep wiU be similar to that of HIV.16|) . 
First, we make the following substitution 

M{p)=e^P>'M{p), -l<^i<l, (IV.33) 



(IV.29) 
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in order to eliminate the p^/i^ term in (jIV.17|) . We obtain 



[(1 - /i2)M']' ± 2p(l - /i2)M' 



(-2i?2 T 2p)/. + - A) 



M = 0. 



(IV.34) 



To simplify notation, let us just consider the case M(/i) = e P^Af(/i), but note that for Af = eP'^A/(/i), we need only 
make the changes of p — > — p in IIIV.37|) below. Write 



(IV.35) 



fe=0 



where P^{fi) are the associated Legendre polynomials, and substitute (|IV.35|I into (|IV.17|) . We obtain a 3-term 
recurrence relation 



Ctnfn-l - Pnfn + 7"/"+! = 0; 



n = 0,l,2,...;/_ 



0, 



(IV.36) 



where 



1 



2(m + n)-l 
/3„ = A — + (m + n){m + n + 1) 
2m + n + 1 



In 



2{m + n) + 3 



2ni?2 + 2pn{m + n)], 
^)(m + n + 1), 
{-2i?2 - 2p(m + n+ 1)}, 



(IV.37) 



and, consequently, again the continued fractions of the same form as (jIV.29|) . The continued fractions obtained here 
should be coupled with the continued fraction ljIV.29p for the variable /i to solve A and p. 
In the homonuclear case, i?2 = R{Za — Z}j)/2 = 0, equation IIIV.17P reduces to 



[(i-^2)^fT+ -^ + pV- 



In this case, several different optional representations of M can be used: 



(a) M{fi)^{l^fi^)'^ Eckf^' 



M = 0. 



k=0 
oo 



k=0 



(b) A^(m) = (1 - m') V ^ CfeP™+2,(A.), M{p) = (1 - A*^) V ^ CfcP™^2,^i(M); 



fc=0 



k=0 



(IV.38) 
(IV.39) 



(c) Af(/x) = e±P'^(l - E T m)' 

/c=0 



In Appendix ^ we discuss expansions of solution near A ~ 1 and A ^ 1 and their connection with the James- 
Coolidge trial wave functions. 

As a conclusion of this section, we note that the eigenstates of the hydrogen atom given in the preceding subsection 
can also be easily represented in terms of the prolate spheroidal coordinates. We let the nucleus of H (i.e., a proton) sit 
at location a where (0, 0, —R/2) with Za = 1 while at location b where (0, 0, R/2) we let Zb = 0. Thus, the hydrogen 
atom satisfies Eq. (jIV.13|) in the form 



Hip ■ 



Now, in terms of the prolate spheroidal coordinates (|A.3(I in Appendix A, and 

ip{X,fi,(l)) = A{X)M{fi)<^{(f>), where <^{(f>) = e™"^ 
in the form of separated variables, we have 



(IV.40) 



(IV.41) 



1 



d_ 

2 1 



Af - 



d_ 



A 



(A2 - /i2)^2 

(A2-l)(l-/i2) 



AfA 



RX- IJL 



MA = EM A, 
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which has two fewer terms than ljA.6|) does as now Z\, = 0. Set = —R^E/2, we again have (|IV.16p and ljIV.17|) 
except that now Ri = R2 = R/2 therein. The rest of the procedures follows in the same way with some minor 
adjustments as noted above. 

The above discussion also leads to a sequence of identities between (|IV.6ll and l|IV.4ip . as 

V'^^^ {x, y, z) = t/i^^^ y^z+-^ ' 

where tp'^^'^ (x, y, z) is an eigenstate of the hydrogen atom obtained from (|IV.2|I but expressed in terms of the Cartesian 
coordinates {x,y,z) while ijj'^-^'' {x , y , z) is that for the solution of HIV.40p . 



C. The many-centered, one-electron problem 



When iVe = 1 and > 3 in (|1V.1|I . we have a molecular ion with three or more nuclei sharing one electron. 
A simple example is a C02-like structure, with A^ = 3. For such a problem, separable closed-form solutions are 
extremely difficult to come by from the traditional line of attack. However, we want to describe an elegant analysis 
by T. Shibuya and C.E. Wulfman [s^ (see also the book by B.R. Judd |63|) which works in momentum space and 
expand electron's eigenfunction as a linear combination of 4-dimensional spherical harmonics. This analysis may offer 
useful help to the modeling and computation of complex molecules after proper numerical realization. 

The model equation reads 



N 

2' ^,\r-RA 



4v'-y^^^ U(0 =^VXr), r={x,y,z)GK', (IV.42) 



where Rj are positions of the nuclei. 

Appendix IeI shows how to reduce the problem to a matrix form. Here we provide the answer for the energy; it is 
determined from the solution of the eigenvalue equation 

Pc = V^2Ec, (IV.43) 

where c is an infinite dimensional vector, P is an infinite matrix with entries 

[III 1*1 
S'n^i^m" (R-i ) - Sn"7''m" i^j ) J (IV.44) 



n, n\ n" = 1,2,3,...;^, £' , £" = 0, 1, 2, . . . ; m, m', to" = + 1, . . . , -1, 0, 1, . . . , £ - 1, 

the matrix S is given by an integral over 4-dimensional unit hypersphere S3 with the surface element dfl 
sin^ X sin 6dxd9d(j), 



SllfP^iiRj)^ [ expiiRj ■p)Ynimifl)Yn'l',n'ifl)dn, (IV.45) 



'S3 

Ynemi^) is a product of the spherical function Yira{0, 4>) and the associated Gegenbauer function C^(x) 
The 3-dimensional vector p in Eq. (jIV.45ll has components 



= P sin cos 0, pj, = p sin sin 0, pz—pcos9, where p — V— 2i?tan(x/2). 

In practice, the infinite matrix P in ljIV.44p is truncated to a finite size square matrix according to the quantum 
numbers {nim) for which the restriction n < np is specified for some positive integer no- 

In the derivation, if we restrict N = 1, Zi = 1 and set J?i = 0, then the matrix P is diagonal and we recover the 
hydrogen atom as derived in Subsection llV Al Obviously, if A^ = 2, by setting Ri = (0, 0, — i?/2) and R2 = (0, 0, R/2), 
we should also be able to recover those H^-like solutions given in Subsection llVI B. 

The 1-electron one-centered or two-centered orbitals derived in Subsection llVI A and II VI B will be utilized frequently 
in the rest of the paper. At the present time, there is very limited knowledge about the 1-electron many-centered 
orbitals as discussed in Subsection llVI C . There seems to be abundant space for their exploitation in molecular modeling 
and computation in the future. 
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V. TWO ELECTRON MOLECULES: CUSP CONDITIONS AND CORRELATION FUNCTIONS 



A. The cusp conditions 



In the study of any linear partial differential equations with singular coefRcients, it is well known to the theorists 
that solutions will have important peculiar behavior at and near the locations of the singularities. We have first 
encountered such singularities in Subsection llll C"T1 Here we give singularities of the Coulomb type a more systematic 
treatment. The critical mathematical analysis was first made by Kato |64l | in the form of cusp conditions for the 
Born-Oppenheimer separation. 

Consider the following slightly more general form of the Schrodinger equation for a 2-particle system 



where 



{H-E)^ = Q, (V.l) 

= V?-— V2-^-^-^-^ + ^ + M!i. (V.2) 

2toi 2m2 Via rib ''2a ri2 R 

The operator H has five sets of singularities, at 

ria = 0, rib = 0, r2a = 0, r2b = and ri2 = 0. (V.3) 

It has been proved by Kato "g^ that the wave function ^ is Holder continuous, with bounded first order partial 
derivatives. However, these first order partial derivatives dtp/dxi, etc., i = 1, 2, . . . , 6, are discontinuous at ljV.3|) . In 
the terminology of the mathematical theory of partial differential equations, ljV.2(l is said to have a nontrivial solution 
in the Sobolev space iJi(R'5). 

We now discuss the cusp conditions at these singularities. What is a cusp condition! It can be simply explained in 
the following paragraph. Let us elucidate it for the two particle Hamiltonian (|V.2p : for a multi-particle Hamiltonian 
the idea is the same. 

In order for the wave function ip to satisfy the eigenvalue problem (|V.1|I at the singularities (|V.3|l , the kinetic energy 
operators — V^/2toi and — V|/2to27 after acting on ip, must produce terms that exactly cancel those singularity terms 
in the potential in order to give us back just a constant E times "0, because the wave function ip is bounded everywhere 
in space, including the points where the nuclei are located, without exception. One can see that, if the cusp conditions 
are not satisfied, then there is some unboundedness at the singularities l|V.3|) which can affect the accuracy in numerical 
computations. Conversely, if the cusp conditions are satisfied, this normally improves the numerical accuracy. 

In case we don't know the exact eigenstate, but only a certain trial wave function, say (f>, then {H — E)(/) = will 
not be satisfied in general. Rather, we have 

{H - E)(t)(ri,r2) = /(ri,r2) 

for some function / depending on the spatial variables ri and r2- However, we can insist on choosing parameters 
in (f) such that the residual /(ri,r2) is a bounded function everywhere] in particular, /(j"i,r2) cannot contain any 
singularity at ljV.3p . We say that the trial wave function (j) satisfies 

(i) the electron-nucleus cusp condition at a (resp. b) if / is not singular at a (resp. b); 

(ii) the interelectronic (or electron- electron) cusp condition if / is not singular when ri2 = 0. 
For example, in the simple case of a hydrogen atom. 



let 4>{r) = Ce""'' be a trial wave function. Then for any E, 



1 



{H-E)(b^Ci ^^"^ °' -iE + o?l2)e- 



The singularity l/r can be eliminated only by choosing a = I. This is the cusp condition, which actually forces to 
be the ground state (with E ~ —1/2). The profile of 0, as shown in Fig. 1101 illustrates the appearance of a cusp at 
the origin. 
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FIG. 10: A 1-dimensional cross section of 4>{r) = e *" showing a cusp. 



In Appendix IFl we derive the cusp conditions for the two particle electron wave function ip of (|V.ip : 



dip 



dria 



-miZaip{ria = 0), 



ria=0 



dip 



dru 



-miZbip{rib = 0), 



(V.4) 



dr2a 



-m2Za1p{r2a = 0), 



dip 



dr- 



'2b 



-m2Zbip{r2b = 0), 



1-26 = 



(V.5) 



dip 



dr\2 



mim2 



qiq2'ip{ri2 = 0). 



(V.6) 



Eqs. ljV.4|) . ljV.5|) are the electron-nucleus cusp condition, while Eq. (|V.6|1 is the interelectronic condition. 

In forming trial wave functions from one-centered or two-centered orbitals for a homonuclear diatomic molecule a 
commonly used wave function is 



'>i^{ri,r2) 



{ri)cP{r2)f{ru), /(ria) = 1 + -ri2, 



(cf. gJ7\ 



(V.7) 



where (p{ri), i — 1,2, is an orbital for the molecular ion. In Appendix IFlwe show that l|V.7|l satisfies the interelectronic 
cusp condition. If, however, (pi ^ (p2, then the trial wave function 

1 



V'(»"i,»"2) = <Piiri)(p2ir2) (^1 + -^^2 

satisfies the interelectronic cusp condition if and only if 

Mr)VMr) - 0i(r)V02(r) = 0. 

The actual verifications of cusp conditions for specifically given examples of trial wave functions in the cases of 
one-centered orbitals or their products are not difficult. But such work is nontrivial when the trial wave functions are 
expressed in terms of prolate spheroidal coordinates. In Appendix |3 we illustrate through concrete examples how to 
carry out this task. 



B. Various forms of the correlation function f(r\2) 



We have learned the importance of the interelectronic cusp condition in Section FV Al But there are, in addition, 
three important constructs that are crucial for diatomic calculations: orbitals, configurations and electronic correlation. 
In this section, we compile a list of often cited correlation functions / which help the satisfaction of the interelectronic 
cusp conditions in the context of Eq. (|V.6|) . The study of the other two, i.e., orbitals and configurations, will be 
addressed in the next section. 
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(1) /(ri2) = l + iri2. 

This is the simplest possible interelectronic configuration function. The specific form is due to the correlation cusp 
condition only. We derive it as follows (see Patil, Tang & Toennies Consider two charged particles which are 

described by the Schrodinger equation 

^V? - -^V^ + ^) ^ = e^. (V.8) 

2toi 2m2 ri2 / 

First, transform the above equation to the center-of-mass coordinates (cf. Appendix IH)! 

V 2M 

where 



nr , mim2 miri+TO2r2 n^1n^ 

M = mi +7712, fi= , S= , ri2=ri-r2. (V.IO) 

mi + m2 m-i + TO2 

In the ground state ip is independent of S. Near the singularity point ri2 = 0, the wave function -ip has a local 
representation as a power series 

i, = Co + C,n2+0{rl^). (V.ll) 
Substituting IjV.lip into (|V.9p and using = ^^^^ (''^ J^) obtain 

— (-— + qiq2Co] + (-— + Cmq2) +--- = e{Co + Ciri2 + 0(r?2)) (V.12) 
''12 V M J \ fJ- ) 

The above mandates that the coefficient of l/ri2 must vanish, that is 

Ci = i-iqiq2Co, 

which for fi = 1/2 (i.e., mi — m2), and qi = q2 = I, yields 



This gives us the small ri2 behavior 



Ci - ico. 



^ = Co[l + \r^2 ) , (V.13) 



where we have dropped all the 0{r\2) terms. The asymptotic expression (jV.13(l motivates the choice of f{ri2) in such 
particular form. 

This simple f{ri2) captures short distance interelectronic interaction very well. It offers elegant representations of 
molecular orbitals and great facility to computation. Nevertheless, its asymptotic behavior of linear growth for large 
ri2 is not physically correct. 

When we write molecular orbitals as 



'P{ri,r2,ri2) = (?!)(ri,r2)/(ri2), 

if the function (j){ri,r2) is already quite small in the region where ri2 becomes large compared to 1, then this simple 
f{ri2) = 1 + ^ri2 can work quite well (Kleinekathofer et al. [13, PP- 2841-2842]). 

(2) /(ri2) = l + ^e='--/'', (d>0) 

This function was proposed by Hirschfelder [g^ where d is a variational parameter. Its profiles is shown in Fig. It 
satisfies the cuspcondition near ri2 = 0. This function was used by Siebbles, Marshall and Le Sech [g^. Nevertheless, 
Le Sech et al. |^ |66] reported that in performing variational calculations by writing /(ri2) = 1 + ^e""''^^, they 
found that a is computed to be very close to for small and intermediate R. That is, it virtually degenerates into 
/(ri2) = 1 + ^ri2- However, at large R one should use a ^ in order to obtain the correct dissociation limit. 



40 




FIG. 11: Graph of f{rvi) = 1 + ""i^/"*, where the maximum happens at rn = d. 



lYvoTn physical considerations, there is no reason to believe why /(ri2) should have a local maximum as shown in 
Fig. ^] Even though the choice of this /{ru) seems satisfactory asymptotically for both ri2 small and large, it may 
not be satisfactory for medium values of ri2- 

(3) = 1-3-^6-^'--, (A >0) 

This correlation function is partly motivated by Hirschfelder's work ^6^, and partly by Hylleraas study of the helium 
atom ^SQJ- It was introduced by Kleinekathofer et al. [37i | . 
At small ri2 we have the expansion 



f^r,2) = ( 1 + -ri2 - -rl^ ± ■ ■ 

■'^ ' 1 + 2A V 2 4 1^ 

Therefore the cusp condition is satisfied for any A > 0. This A can be either used as a variational parameter, or be 
determined from a given Hamiltonian. For example, for a helium-like 2-electron atom with nuclear charge Z, using a 
perturbation argument, Kleinekathofer et al. 37] analyzed that the best value for A is 

A^lz-i. 
12 3 

This f{ri2) has a monotone profile and correct asymptotics for both ri2 small and large. See Fig. 1121 



f 




^12 

FIG. 12: Graph of /(n2) = 1 - ir^x e"'^'"'^ • 



(4) /(ri2) = e^'--. 

This f{ri2) satisfies the correlation-wave equation 
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with a negative energy —1/4. For small ri2, its expansion is 

Therefore, it satisfies the correlation-cusp condition and its asymptotics for small ri2 is good. However, this f{ri2) 
has exponential growth for large and is thus physically incorrect. 

(5) ) = !l^.£2(il^ 

tr ri2 

where Fj(ri,p) is the Coulomb wave function regular at the origin, j is an integer, t,k are separation of variables 
constants and r = |ri + r2|. 

This is perhaps the most complex form of the correlation function in the literature, given by Aubert-Frecon and 
Le Sech hj. It comes from solving 



Rewrite the above in center-of-mass coordinates: 



Assume that / depends only on r, ri2 and 7, where 7 is the angle between S and ri2 (see Fig. I13|l . 




A B 

FIG. 13: The variables r, ri2 and angle 7 in the center-of-mass coordinates. 



(V.14) 



(V.15) 



Then equation IjV.lSjl can be written as 



df 



1 d 



dri2 V' ^^'9ri2 



1 d ( 2d.f 



dr \ d' 
1 



ri2 dq oq 



+ —f = ef, 
ri2 



q = cos 7. 



dq dq 
Equation HV.16p can be separated by writing 

/ = P,(gV(rK(ri2), 

where 

Pj is the Legendre polynomial of degree j and g^{r) and u^{r 12) satisfy, respectively, 



I- 



2 d j{j + 1) 
r dr 



d^ 



2 d 



dr{2 ri2 dri2 



j{j + 1) 



7-12 



+ en >u^{r 12) = 0, 



(V.16) 



(V.17) 

(V.18) 
(V.19) 
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with 



Take 



Then 



(V.20) 



Po(g)-l, 5?M = 



sinh(<r) o _ ^o(^, fc^: 



where the function Fo(ri,p) belongs to the class of regular Coulomb wave function FL{'r],p), satisfying 



^ 27? L{L + 1) 



for the two parameters 77 and L. 



Note that t and k in ljV.20|) can be used as variational parameters. For small ri2, 



^0 ( ^'^^12 ) = ^ • ^''12 



l + ^r-i2+0((A;ri2)2) 



and, thus, 



7'12 



1+ 2^12+C^((fc^l2)') 



satisfies the correlation-cusp condition for any given k. 

For large ri2, the Coulomb wave function FL{rj^p) has an asymptotic expansion 



Fl 


= gcosOL + /sin^L, 


Ol 


— p — rjln zp — 




Oh 


= argr(L + l- 


hill), 


f 


OQ 

^^fk, g^ 


00 

Y.9k, 




k=0 


k=0 


fo 


= 1, ffo = 1, 




fk+1 


= akfk - ^'fcSfcj 


gk+1 = ikgk + bkfk, 


flfe 


(2fc+ 1)7? 


L{L + l)~k{k + 


(2/c + 2)p' 


^~ i2k + 2)p 



Thus 



r 0/ ^ r Fo{^,kri2) 
hm uj^Xr 12) = Imi — = 0. 



(6) /(r-12) = iFi (-^,2,-2Zri2^ 

where i-F'i(a; 6; a;) is the confluent hypergeometric function satisfying the differential equation 



x—r^ + (b — x)— aw — 0. 

ax^ ax 



(V.21) 

This f{ri2) was given by Patil, Tang and Toennies for the case when the internuclear separation R is large: 

3 



(V.22) 
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Its derivation can be motivated as follows. Consider 



1^.2 ( Z Z Z Z\ 1 

H = V? Vn \ \ \ + \ 

2 2 V^la T\h f2a f26 / Txi R 



(V.23) 



for an H2-like molecule. Guillemin-Zener-type one-electron wave functions |54| suggest the molecular orbital for large 
R: 



^ = (e-^'-i--^^^" + e-^''i''-^''^")/(ri2). (V.24) 
When R is large, electron 1 is localized around A and electron 2 is localized around B, as shown in Fig. 




We have 



FIG. 14: For large R, electron ei is localized near A, and electron 62 is localized near B. 



Tib > ria, r2a > ^26. 



Because of (jV.25|l . we have, from (jV.24|) 



Substituting (06)1 into = Eip, we obtain 

1 



(V? + V^)/ + Z(Vi/) • (fi, + f2b) + — / = O ( ^ 
2 ri2 \R 



(V.25) 



(V.26) 



(V.27) 



Note that Eq. 1V.27|) contains the effect of cross terms. For large R, the electron- electron correlation is most significant 
when the two electrons are colinear and in between the two nuclei. In this situation, either ri2 is antiparallel to ri — r^, 
or ri2 is parallel to r2 — r?,. 

Now, using the center-of-mass coordinates and dropping all the O (1/-R) terms, we obtain 



92 



dr{2 ri2 dri2 



ri2 



dr 



12 



(V.28) 



The solution to (|V.28p . after setting it into the form of (|V.2ip is 



For small ri2, the expansion is 



/(ri2) = iFi (- — ,2,-2Zri2 



/(»-12) = 1 + 2^12 ^ ^12 + ' 



Therefore, the correlation-cusp condition is satisfied. For large ri2, the asymptotics is 

/(ri2)~C(-2Zri2)* 
(cf. Abramowitz and Stegun [gI p. 508, 13.5.1]). 
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VI. MODELING OF DIATOMIC MOLECULES 



In this section, we give a survey of major existing methods for the numerical modeling of diatomic molecules. These 
methods provide approximations of wave functions either in explicit form through properly selected ansatzs, or in 
implicit form through iterations as numerical solutions of integro-partial differential equations. 

The methods and ansatzs to be described below are 

(1) The Heitler-London method: 

(2) The Hund-MuUiken method; 

(3) The Hartree-Fock (self- consistent) method; 

(4) The James-Coolidge wave Junction; 

(5) Two-centered orbitals; 

Items (1)~(3) above historically are associated with one-centered orbitals, while (4)~(5) are based on two-centered 
orbitals. But this dichotomy is not inflexible. An example is a hybrid type containing both one-centered and two- 
centered orbitals considerd in Subsection IIIII C .6 . We now discuss them in sequential order below. Each approach has 
a set of modeling parameters which can be optimized through calculus of variations. In particular, we will point out 
what these parameters are. 



A. The Heitler— London method 



This method has the longest history. It was developed by Heitler and London during the 1920s soon after Heisenberg 
laid the quantum mechanical foundation of ferromagnetism. The method is usually called the valence-bond (or atomic 
orbital) method. In this method, each molecule is thought of as composed of atoms, and the electronic structure is 
described using atomic orbitals of these atoms. 

Here we present a version of refined Heitler-London approach due to Slater f6§j. In the method, electron spin- 
orbitals are taken from a determinant 



(VI.l) 



which satisfy the Fermionic property of the Pauli exclusion principle. The orbitals in (|VI.1(I are called the Slater 
orbitals. Let us consider the 2-electron case, i.e., n = 2 in (jVI.ip . The spin-orbital Ui{j), i,j = 1,2, consists of 
(i) the electron-atomic orbital part 



U,{1) 


Ul{2) ■ 


• Ui{n 


«2(1) 


"2(2) • 


■ U2{n 


Mn(l) 


U„(2) 


Un{n 



a(l) = \/— e" 

TT 



6(l) = \/-e 

TT 



a(2) = \ — e~ 

TT 



K2) = \/— e 



(VI.2) 
(VLB) 



where in ljVI.2|) and l|VI.3|l . the atomic electron wave functions are centered at, respectively, a and b. 
(ii) the spin part 



spin a(l), a(2), a{j) ^ \s,ms), t, fori = l,2, 
spin /3(1), /3(2), f3{j)^\s,m,) i, fori = l,2. 
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The linear combinations of the total spin-orbital wave function that are antisymmetric are tabulated below: 





spin-orbital wave functions 


Ms (total spin) 


singlet state 


[a(l)6(2) + 6(l)a(2)][a(l)/3(2) - /3(l)a(2)] 







Kl)6(2)-fe(l)a(2)][a(l)a(2)] 


1 


triplet states 


[a{l)b{2)~b{l)a{2)][a{m2)+(3{l)a{2)] 







[a{l)b{2)~b{l)a{2)mi)Pm 


-1 



(VI.4) 



Table 1. Spin-orbital wave functions of singlet and triplet states 

The singlet state has lower energy than the triplets. For example, if we aim to calculate the ground state of H2, we 
use ljVI.4|l as the trial wave function to minimize the total energy 

subject to (*|*) = 1, 

where 



1^2 1^,2 1 1 1 1 11 
H = — Vi V2 \ h — 

2 2 ria ru r2a 7'26 ''12 R 

^ ^ M[a{l)b{2) + b{l)a(2)], TV = a normalization factor. 



(VI.5) 

(VI.6) 
(VI. 7) 



Heisenberg, and Heitler-London's basic point of view is that H in IIVI.6I) can be written as 



* 2 1 ria 



1 

?'2b 



1111 

rib r2a ri2 R 



-7:^1- — ] + I-- - + — + 



where the terms inside the third pair of parentheses above can be viewed as a "perturbation." 

Since the Heitler-London method is quite fundamental in molecular chemistry, let us give some details about the 
calculation of ljVI.5|) given l|VI.6|) and ljVI.7|) . 

Define 



S = the overlap = (a(l)|6(l)) 



~aria,„-arib 



dx ■ 



{dx — dxidx2dx'i) 



-olR 



1 + aR 



Then (a(l)6(2)|a(2)6(l)) = (a(l)|6(l))(6(2)|a(2)) = S^, and the normalized state for the singlet or the triplets, without 
spin, is 



* = 



a(l)fo(2) ±a(2)&(l) 



such that (*|*) = 1. 



v/2(l±52) 

The total energy of the singlet and triplet states can now be written as 



E± 



^ 'a(l)6(2)±a(2)&(l)|-iv?-iv^ ^ 



1 1 



2(1 ±52) 



ria rib r2a f2h ri2 R 



a(l)6(2)±a(2)6(l)). (VI.8) 



The integrals involved are given in Appendix ^ Using these integrals, we are able to write down the total energy 
E = KE + PE as follows: 

KE± = kinetic energy 

= 2(i±^2) Hl)^(2) ± a(2)&(l)| - iv? - \vlHm2) ± a(2)6(l)) 



1±S'2 



{1^2KStS'^)] 
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PE± = potential energy 



^ (a(l)6(2) ± a(2)&(l)| ~- - + — + 4|a(l)K2) ± a{2)b{l)) 



2(1 ± ria rib »'2a r2b ri2 R 



a 



(1±52) 



(-2 + 2 J + J' ± 4i4:5' ± iT') + - 



where w = aR, the other symbols are defined in Appendix The parameter a can be used as the variational 
parameter to minimize (|VI.5|I (691] . 

Consider the special case when a = 1. Then 



E± = -l 

where 



Ho = [[ a'(l)6'(l) (-- - + — + ^] dxdy 

JJ V ^Ib f2a ri2 R; 

= 2 J + J' + 4 

R 



is the Coulomb integral, and 



Hi= I [ a(\)h{\)a{2)b{2)(-- L + ± + 1^ dxdy 

JJ V ''lb '"la ''12 -K/ 

c2 

= 2KS + K' + — 
R 

is the exchange integral. According to numerical values computed in Slater '6^, Table 3.2], e.g., it is known that Hi is 
usually many times larger than Hq, and is largely responsible for the attraction between atoms in forming a molecule. 

In Fig. El we plot the ground and first excited '^S+ state potential energy curves E{R) of the H2 molecule. 
When a = 1 the ground state curve yields the binding energy of 0.116 a.u.=3.16 eV; the value must be compared 
with 4.748 eV obtained by Kolos and Roothaan When a (effective charge) is treated as a variational parameter 
[70| the calculation yields the binding energy of 0.139 a.u.=3.78 eV and the bond length of 1.41 Bohr radii. For the 
■^S+ state the effective charge and the a = 1 curves are practically indistinguishable. 



B. The Hund-Mulliken method 



In 1927, Robert MuUiken worked with Friedrich Hund and developed the Hund-Mulliken molecular orbital theory 
in which electrons are assigned to states over an entire molecule. Hund-Mulliken 's molecular orbital method was 
more flexible and applicable than the traditional Valence-Bond theory that had previously prevailed. Because of this, 
MuUiken received the Nobel Prize in Chemistry in 1966. 

The approach has some similarity to the Heitler-London's, so we can inherit the notation from there. Its special 
feature is that a linear combination of the molecular gerade (g) and ungerade (u) states are used: 



ly^ 



3v- 



(9+9-) ■■ 

(u+u^) : 

■■ {9+U+) 
'u- +g-'> 



a(l)+fc(l) a(2)+b(2) Tli2-ilT2 

V2(TTs) V2(TTs) 72 ■■ 

a(l)-fc(l) a(2)-b(2) Tli2-ilT2 
V2(T^ V2(T3S) V2 



/2(1-S2) 
a(l)b(2)-b(l)a(2) 
2(1-S2) 



V2(l-S') 
a(l)b(2)-b(l)a(2) 

V2(l-5") 



%/2 ' 






T1T2, 


(M, = 


1) 


2+i.lT2 

V2 ' 


{Ms = 


0) 


ili2, 


(Ms = 


-1 



(VI.9) 
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FIG. 15: Ground and first excited state energy E{R) of H2 molecule for the Heitler-London wave function (solid lines) and 
the "exact" energy of Ref. P| (dots). 



with 



a{i) = 



aii) + hii) aii) ~ bii) 

g: , , u: 



z = l,2, 
i = 1,2. 



Especially, note that the first three states in ljVI.9p . i.e., ^Sg- and ^Sg- signify the possibihty of double 

occupancy of the two electrons at a single nucleus as products of a(l)a(2) and 6(1)6(2) appear in the wave functions. 
Such states, chemically, represent ionic bonds. On the other hand, the last three states in l|VI.9|l i.e., '^S^, '^E~, 
agree with the triplet states in (|VI.4|I . 

We denote the six molecular orbitals in IjVI.Qp in sequential order as ^ti, ^"2, . . . , ^e- Then the energy of any linear 
combination E^^j^Cj^Pj corresponds to a quadratic form: 

/ 6 6 \ 6 

\j=l k=l I j,k=l 

where Hjk are the (j, fc)-entry of the following symmetric matrix 





H12 


H12 


H22 



o 



Hi 



O 



3E„ 



H. 



(VI. 10) 



Uii = 



n 



H22 = \ s 



n 
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Specifically, 



1 =p 5 T 2if 



libs' 



Hi 



-2 + 2J±4ii' 
lis* 



J' 



J' + 2X' ± 4L 1 



2(1 ±5)2 



2(1-^2) ; ' 



n 



1 + 2KS + 



1-^2 
1 + 2KS + 52 



-2 + 2J-4i^5 J-2K' 



l-S'2 



1 

l-S'2 ' 2(1 -S2) ' W 

-2 + 2J + J' -AKS -K' , 1 

w 



(1-S2) 



where 



L = [l/a) / a2(l)a(2)6(2) — dxdy = e" 

ri2 



w -\ e 

8 IQw 



-2w 



The other symbols are defined in Appendix [J 

In the calculation of the ground state energy, the sub-block of 2 x 2 matrix in the upper left corner of ljVI.10|) 
plays the exclusive role as the remaining diagonal 4x4 block in ljVI.10|) contributes no effect. We thus determine the 
ground state energy E by the determinant 



Hii — E H12 
H12 H22 — E 



0, 



E^ 



Hii + H22 ± \/ {Hii — H: 



The value E- will correspond to the ground state energy. 

For the Hund-Mulliken method discussed above, again a is the variational parameter. In Fig. ^| we plot the 
ground and first excited state potential energy curves E{R) of the H2 molecule. When a = 1 the ground 
state curve yields the binding energy of 0.119 a.u.=3.23 eV. The effective charge calculation yields the binding energy 
of 0.148 a.u.=4.03 eV and the bond length of 1.43 Bohr radii. The curve is identical to the Heitler-London E{R) 
(see Fig. [T3|) . 



C. The Hartree— Fock self-consistent method 



This is perhaps the best known method in molecular quantum chemistry and it works for multi-electron and multi- 
center cases. In computational physics, the Hartree-Fock calculation scheme is a self-consistent iterative procedure 
to calculate the optimal single-particle determinant solution to the time-independent Schrodinger equation. As a 
consequence to this, while it calculates the exchange energy exactly, it does not calculate the effect of electron 
correlation at all. The name is for Douglas Hartree, who devised the self consistent field method, and Vladimir Fock 
who reformulated it into the matrix form used today and introduced the exchange energy. 

The starting point for the Hartree-Fock method is a set of approximate orbitals. For an atomic calculation, these 
are typically hydrogenic orbitals. For a molecular calculation, the initial approximate wave functions are typically a 
linear combination of atomic orbitals. This gives a collections of one electron orbitals, which due to the Fermionic 
nature of electrons must be anti-symmetric; the antisymmetry is achieved through the use of a Slater determinant. 

Once an initial wave function is constructed, an electron is selected. The effect of all the other electrons is summed 
up, and used to generate a potential. This is why the procedure is sometimes called a mean-field procedure. This gives 
a single electron in a defined potential, for which the Schrodinger equation can be solved, giving a slightly different 
wave function for that electron. This process is then repeated for all the other electrons, which complete one iteration 
of the procedure. The whole procedure is then repeated until the self-consistent solution is obtained. 

Here we discuss the method in detail. Consider the Hamiltonian of a multi-electron and multi-center molecular 
system (with n electrons and N nuclei, respectively) under the Born-Oppenheimer separation in the following form 

j = l l<j<k<n ■' j = l k=l ' ' 
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FIG. 16: Ground and first excited state E{R) of H2 molecule for the Hund-MuUiken wave function (solid lines) and the "exact" 
energy of Ref. ;9| (dots). 



For a closed-shell system, all electrons are paired up and n must be an even number. The trial wave function is chosen 
in the form of a Slater determinant tp = \\(f>i{ri)4)i{r2)(f>2{i'3)(f>2if4)---4'n/2{'''n-i)4>n/2i'''n)\\, where || . . . || denotes 
determinant and each (j>i corresponding to a Ui in (6.1) is a molecular orbital (MO) and the one without a "bar" 
on top denotes a spin-orbital with a-spin (spin- up), and the one with a "bar" on top denotes a /?-spin (spin-down) 
orbital. We assume that the MOs (the spatial part) are orthonormal: 

/ (t'*i{r)(f)j{r)dr^Sij, 

and we stipulate this orthonormality in our subsequent calculations. For the ground state calculation, we use the Ritz 
variational method. Substituting the trial wave function into {ip\H\'ip)^ we have 

n/2 

{m\^) = 2Y,hu + 5] 5](2J„- - /C„), (VI. 12) 

i=l i j 

where, 

N 



Minimizing (jVI.12p subject to = 1 leads to a set of n/2 Hartree-Fock equations for the spatial molecular orbital 

0,(r): 

Fcj,,{r) =\h + ^[2J, - Kr] I 0,(r) = £,</),(r), j = 1,2,..., n/2, (VI.13) 
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with h being the single electron Hamiltonian 

N 



Jr being the Coulomb interaction operator defined by 



and Kr being the energy exchange operator defined by 



KrMr) - |y^^ Kir')^jirf)j-—^^drf'i Mr)- 

F, called the Fock operator, involves those unknown MOs ((/)j(r)'s) complicatedly (note that summations in F are 
over all occupied orbitals, including (j>j{r) itself); this is not a simple eigenvalue problem. A self-consistent field (SCF) 
method is needed to solve the Hartree-Fock equation. For that, we take all MOs in F to be known by guessing a set 
of MOs and plugging them into F. Then, the Hartree-Fock equation becomes a normal eigenvalue problem. We solve 
this problem and compare the resulting MOs (eigenvectors) to those MOs we substitute into F. If they are different, 
we put these new MOs back to F and do the same calculation again. This is done iteratively until all 0j's converge. 

The problem is then how do we guess (or express) those MOs. An often used method is what we called the linear 
combination of atomic orbitals (LCAO) method. In this method, we expand the unknown MOs linearly in a fixed 
basis set, such as 

— ^ C'^"Xai('")i e C; C is the complex number field. (VI. 14) 

Here, Xfj-ir) can be any functions deemed appropriately (not required to be orthonormal), and they are often approx- 
imations to atomic orbitals with respect to the individual centers. We still refer them as atomic orbitals (AOs) in 
the following. By multiplying (p* (r) on equations ljVI.131) and integrating with respect to r, we can rewrite those HF 
equations (|VI.13|I as 

F.j^eAj, (VI. 15) 

with 

F,, = h,, +Y,['2i^j\rr) ~ Hjr)]. 

r 

If we then substitute expansions ljVI.14|l back into the new HF equations ljVI.15|l , we obtain 

J2 ^A.-^- = J2 St..C,,e,, (VI.16) 

where 



The whole procedure for the Hartree-Fock Self-consistent Field (HF-SCF) calculation can be summarized as follows: 

(1) Choose a basis set {x^}- 

(2) Calculate integrals over this basis set (and store in memory). 

(3) Guess the MO coefficients C^i. 
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(4) Construct D and then F. 

(5) Solve equations (jVI.16|l for new C^;, and iterate until convergence results. 

Note here the total electronic energy can be calculated by use of equation IIVI.12p . 

Good choice of the basis set is very important to the HF-SCF method. Otherwise, basis set can fail the whole 
calculation, produce bad or wrong results or make calculations become very expensive. A natural choice of the basis 
are hydrogenic wave functions, or Slater-type orbitals (STOs) 

xf™ = exp{-'CrAyAYi,n{OA,<l>A), 

where Yim are the spherical harmonics and the STO is centered on the nucleus A. This basis set works nicely for 
atoms as well as for diatomic molecules (some numerical integrations are required) and linear polyatomics (a large 
number of numerical integrations need to be done). No algorithm for nonlinear molecules has been developed till now. 
Another basis set is the Gaussian-type orbitals (GTOs). This basis set is proposed by S.F. Boys [tJ in the following 
form 

x[;™ = exp{-ar\)r\Yir„{eA,(l}A), 

or in Cartesian form, 

This basis set has been widely used since those multicenter integrals in F can be easily done fjl| over this basis set. 
In principle, it works for all molecules. However, there is a trade-off. GTOs are totally conjectural functions which 
have no concrete physical meanings. Because of this, GTOs are seldomly used to form basis sets directly. They are 
often used to estimate STOs (by expressing one STO into linear combination of several GTOs) and then form the 
basis set, such as the ST0-3G (3 GTOs are used to estimate one STO) basis set. 

So far we have discussed the HF-SCF method for closed-shell systems. For the open-shell systems, there are 
unpaired electrons, orbitals could be doubly occupied or singly occupied. This makes calculations more complicated. 
The simplest way to deal with this problem is to treat one doubly occupied orbital as two independent orbitals. For 
example, orbitals (f>i and 0i are treated independently and they are not required to have the same spatial part. With 
a few modifications, the closed-shell HF-SCF method can be migrated to this case, which is called the unrestricted 
open-shell HF-SCF method. If we do put restrictions on and requiring them to have the same spatial part, then 
we need to treat the doubly occupied and the singly occupied orbitals differently. Calculations can be done still in a 
more complicated way called the restricted open-shell HF-SCF method ,72j. 

Many strategies have been developed to improve the HF-SCF calculations such as introducing unoccupied molecular 
orbitals (virtual orbitals), or using two (double-zeta basis set) or three (triple-zeta basis set) STOs to describe one 
molecular orbital (note here that STOs could be linear combinations of GTOs). 

The HF-SCF method succeeds in many purposes of calculations such as the ground state energy, chemical bond 
lengths, angles calculations, etc. It yields 3.636 eV for the binding energy of H2 molecule 9] (to be compared against 
the "exact value of 4.745 eV). Nevertheless, this method fails to describe the long distance behavior of a chemical bond 
(see Fig. I17|l . More subtly, the HF method views that each electron interacts with a mean potential field generated 
by the other electrons. It does not take into account the electron correlation. Researchers have been trying to mend 
these by combining other methods, such as the multiconfigurational SCF, configuration interaction and interelectronic 
correlation terms in their calculations. 



D. The James— Coolidge wave functions 



James and Coolidge |20| suggested a wave function of the form 

00 

^ = —e'^^^^+^^^ E c^n,kp{XT>^2^^i^^2 + >^l^2^^l^^i)p' (vi.i7) 

= 1 

for the H2-molecule, where ip is a. function of five variables 
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FIG. 17: Ground state E{R) of H2 molecule calculated by the self-consistent Hartree-Fock method (solid line) and the "exact" 
energy (dots). 



^From homonuclear (such as H2) symmetry considerations, j + fc in (|VI.17|I must be even. But this restriction can be 
removed for heteronuclear cases. For given i? > 0, the coefficients a, and Cmnjkp for finitely many indices m,n,j,k 
and p, can be used as variational parameters to minimize the energy (iplHl^), subject to the normahzation condition 



{tpltp) = 1. (VI.19) 

For example, for fixed a and R in (IVI.17|) and (|VI.18p . introduce a Lagrange multiplier A for the constraint l|VI.19|) 
and choose only a total of s terms in (|VL17|I by truncation, then the variational problem 

min[(i^|7?|7/^) + X{{^\^) - 1)] (VI.20) 
leads to a set of s linear equations ( 20], p. 826]) for the coefficients Cj, j — 1,2, . . . , s, 

{H12 

< 

. (His 

where 

Hij = {(f>i\H\(l)j), Sij = {(f>i\(f>j), 

(j>i is a typical summand term in ljVI.17() without the coefficient Cj , 
for i,j = 1,2,.. .,s. 

Solving (|VI.21|I then leads to the ground state ipo and its energy Eq. 
The trial wave function (|VI.17|I may be further adapted to 

V' = ^e-"i^^-"^^^e-'^i^i-'3^^^ £ annJkp{XT^2^^{^^2±^l^2^^l^^i)p'' (vi.22) 

= 1 

for the calculation of heteronuclear cases and excited states. 

Next, we address the analytic treatment and provide a complete compendium for the integrals given in (|VI.17|I - 
(|VI.22|I . which constitutes the keystone in the two-centered variational treatment. 



A5ii)Ci + {H12 — XSi2)C2 

XSi2)Ci -t- {H22 ~ XS22)C2 



{His — XSis)Cs — 0, 

{H2s - XS2s)C2 = 0, 



(VI.21) 



A5is)Ci + (i?2s — XS2s)C2 



{Hss — XSss)Cs — 
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The model Hamiltonian that we will be using here is the (homonuclear case) H2 as given in ljVI.6|l . Let us rewrite 
the James-Coolidge wave function as 

*(1,2) = ^a$r(l,2), (VI.23) 

r 

where 

$.(1,2) = ^e-"(^^+^^)Ar'-ArMf M2'-^i2. (VI.24) 

We remark that for the heteronuclear case and for excited states, the trial wave function (|VI.23|) - ljVI.24p can be easily 
re-adjusted to the form ljVI.22p with relative ease. 

In the evaluation of the energy for given two-electron trial wave functions, a typical term is of the form 

Z''{m,n,j,k;e) = i / • ■ • / e-^f^i+^^^A^'A^/xj/i^/ia^^" cos''(0i - <j)2) dXidX2diiidn2d(pid(p2 (VI.25) 



47r2 „ 
where 

M = [(A? - ml 1)(1 - ^?)(1 - (VI.26) 

and v, TO, n, j, k and £ are integers {£ > —1 and others are nonnegative integers). 

Using the above function, we can write every integration in terms of Z'^{m,n,j,k;£). For example, the Coulomb 
interaction energy between the nuclei and electrons may be easily written as follows: 

(VI. 27) 
(VI.28) 
(VI.29) 
(VI.30) 



1 


2 




rai 


i?(Ai + 


Ml) 


1 


2 




ni 


i?(Ai- 


Ml) 


1 


2 






i?(A2 + 


M2) 


1 


2 





rb2 i?(A2 - M2) 

The evaluation of the Coulomb interactions can be done in terms of Z. The evaluation of the Laplacian matrix 
element, however, somewhat more involved. But it also can be expressed via Z'^ . We provide details of the derivation 
in Appendix U 

Since we can write every term in the energy in terms of Z"^, let us express Z'^ via simpler functions defined by 
recurrence relations. Such relations are given in Appendix IK1 For £ > 1 and i/ = 0, one can reduce ^ to or —1 using 
the following identity: 

r\2 = — (A? + A^ + ^2 _^ ^2 _ 2 - 2AiA2AiiM2 - 2Mcos((/)i - ^2)). (VI.31) 
For J/ = and £ = 0, the integration can be evaluated as follows. Given 

Z°(TO,n,j,A:;0) = Z(TO,n,j,/c;0) = ^ /"■/ ^"^"'^'^■^'■'^i'^2 ^^^2 '^Ai dAa d^i ^^M2 rf'/>i ^^2, (VI.32) 

note that the integrals are trivial to evaluate and the /i integrals survive only for even integers j and fc, i.e., 

^ . 2 

[i^d[i = -, for even j ; (VI. 33) 



we arrive at 



A(m:a)A(n:a) , ■ , .in , for even 7 and k, 
Z{m,n,j,k;Q) = I ^ ' ^ ^ ' ' [i3+i){k+i) \ ^ J ' (VI.34) 

1^0, for odd i or k, 
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where we have introduced the function 



A{m 



(VI.35) 



for the A integration. Using integration by parts one can show that the A{m; a) satisfies the recursion relation 



with 



A{m; a) = — [e " + mA{m — 1; a)] , 



A(0;a) = — . 

a 



(VI.36) 



(VI.37) 



Thus we have given a recipe for the evaluation of Z{m, n, j, k; 0) for arbitrary values of the power parameters m, n, j, 
and k. 

Next, consider the case where ^ = — 1, i.e., 

Z(m,n,i,fc;-1) = At / e-2«(^i+^=)A5"A>>J— dAi dAa rf/Ui <^/«2 #i #2 • (VI.38) 

47r^ J J ri2 

Recalling the Neumann expansion for l/ri2: 



1 r, OO T 



r=0 L'=—T 



(r-IH)!' 



(r+H)!, 



p;(A<)g^(A>)p;(/xi)p;(/x2)e 



iv{(t>1—<j>2) 



(VI.39) 



where A< = min{Xi, A2), A> = ma.E(Ai, A2) and P^, Q'^ arc the associated Legcndrc functions of the 1st and 2nd kind 
respectively. After the angular integration, only terms corresponding to = survive, as long as the wave function 
has no angular or ri2 dependence. Separating the A and /x integrals, we arrive at 



Z{m, n,j, fc; -1) = - ^(2t + l)Rr{j)Rr{k)H^{m, n; a) , 
where Rt and Ht are defined by 

Rr{j) = j y Prill) dfl, 

j e-"(^i+^^)A™A^Pr(A<)Q^(A>) d\id\2. 



(VI.40) 

(VI.41) 
(VI.42) 



In the discussion to follow we give recursion relations for the evaluation of the various auxiliary functions. For r = 0, 



Holm, n; a) = A(m; a)F{n; a) + A{n: a)F{ni: a) — T{m^ n; a) — T(n, to; a). 
Here, F{m\a) can be evaluated for arbitrary to by noting its recursion relation 



(VI.43) 



and the initial values 



and 



e-"^A™go(A) dX 



F{m - 2; a) + - [toF(to - 1; a) - (to - 2)F(to - 3; a) - A{m - 2; a)] 



F(0;a) = ^ (In 2a + 7)^ Ei[-2a] — 

a a 



F{l;a) 



(ln2a + 7)e-" ( - + ^ ) - Ei[-2a]e" ( -- + ^ 



a 



a a' 



(VI.44) 



(VI.45) 
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where 



f°° e"* 

Ei(-x) = - / —dt (VI.46) 



and 7 = 0.577216 • • • is the Euler constant; see Appendix I. 

Similarly the quantity T{m, n; a) can be determined for arbitrary values of m, n through 

ml \ ^ a 



with the initial value 



= — \mTini — 1, n; ct) + _F(m + n; 2a)] 
a 

T(0, n; a) = -F{n; 2a) . (VI.47) 

a 

Note that we have so far considered only a special case where r = 0. We turn to the case where r = 1. Once again 
we note the recursion relation for Hilm, n; a): 

Hi{m, n; a) — Ho{m + 1, n + 1; a) — S{m, n + 1; a) — S{n, m + 1; a) (VI. 48) 

where S{m, n; a) can be determined according to 

ml \ ^ a 



1 



with initial value 



mSim — 1, n; a) + Aim + n; 2a)] 

a ' 

S'(0, n; a) -A{n; 2a). (VI.49) 
a 

We now have shown relations needed to evaluate Hr{m, n; a) for particular values of r = 0, 1. Here we summarize 
the evaluation for t > 1 through the following recursion relations: 

Hr{m,n;a) = [(2t - l)'^Hr-i{m + 1, n + 1; a) + (t - l)'^Hr~2{m,n;a) 

- (2r - l)(2r - 3) {Hr-2{m + 2, n; a) + Hr-2{m, n + 2; a)} 

+ 2(2t - 1)(2t - b)Hr-3{m + 1, n + 1; a) 

- (2r - 1)(2t - 7) {Hr-i{m + 2,n; a) + Hr-i{m,n + 2;a)} 

+ 2(2r - 1)(2t - 9)Hr-5{m + 1, n + 1; a) 

+ one of the following: 

{(even t) — (2r — 1) {Ho{m + 2, n; a) + Ho{m, n + 2; a) — S{m + 1, n; a) — S{n + 1, m; a)}] , 
(odd t) + (2t - 1) {2Ho{m + 1, n + 1; a) - S{m, n + 1; a) - S{n, m + 1; a)}] , 

where the starting values Ho{'m, n; a) and Hi{m, n; a) are already given in svrm and (ivr^ . 

For nonzero v, the integral Z'^ can be written in terms of various simple function as discussed below. It is fairly 
straightforward to obtain the recursion relation for the Z'^ function, just by inspection of the equation. For £ > 1, 

Z" (to, n, j, k-J) = Z''{m + 2, n, j, k;£-2) + Z^im, n + 2, j, k;£-2) 

+ Z'' (to, + 2,k;£-2) + Z^im, n,j,k + 2;£-2)- 2Z'' (to, n, j, k;£-2) 

- 2Z'\m + 1, n + 1, j + 1, fc + 1; £ - 2) - 2Z''+^{m, n, j, fc; £ - 2). (VI.50) 

The above recursion relation can be proved in a straightforward manner by using the expansion of iii the elliptical 
coordinates introduced in ljVI.31|) . 
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Higher order terms of are given by 



{m, n, J, fc, - 1) = - 1 £ J^L±^Rl^j)Rl ^k)Hl (m, n; a) (VI.51) 



and 



T^2 



^ oo 

- E(2t + l){Rr{3) - RrU + 2)){Rr{k) -Rr{k + 2)) 



X (Hr{m + 2,n + 2;a) - Hr{m + 2,n;a) 
— Hr{m, n + 2; a) + Hr{m, n; a)). 

For the £ — case, 

Z^(m,n,j,k;0) = 0, 



Z'^{m,n,j,k;0) = {A{m + 2;a) - A{m;a)){A{n + 2;a)- A{n;a)) — 



(j + l)(j + 3)(fc + l)(A: + 3)' 



when j and k are even; otherwise Z^(m, n,j, fc; 0) vanishes. 

The definitions and recurrence relations for i?^ and iJ^ are given below: 



oo /"OO 



i7;(m, n;a)= / e-"(^i+^^) A^(A? - 1)"/2(A2 - 1)^/2f;(A<)Q!;(A>) dAidAa. 



The higher order recursions for H!^{m, n; a) for v — I and = 2 are listed below: 

Hl{m,n]a) = — — — — 7JT-+i(m, n; a) - t{t + l)Hr{m + l,n+ !;«) + — - — ——Hr-i{m,n;a), 
[ZT + L) ZT + l 



and 

(2t+1) 



H^{m,n;a) = -^^^——^H^_^_^{m,n■,a) - (t + 2)(t- l)iJ^(TO + 1, n + 1; a) 



(r + 2)(r + l)2 

+ ^— --j H^_^{m, n; a). 

Zt + L 

The James-Coolidge wave functions may be said to have the most "brawny" power as thousands of coefficients 
Cmnjkp have been calculated with automated computer programs, yielding very accurate values for the binding energy 
of diatomic molecules. Nevertheless, there are certain associated shortcomings: 

(i) The physical insights about molecular bonding seem to be lost; 

(ii) The wave functions in general do not satisfy the correlation cusp condition; 

(iii) For large Ai and A2, the asymptotic conditions are violated (see Appendix Id)| . 
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FIG. 18: Schematic correlation diagram for He, the united atom hmit, and H2, the separated-atom case. 



E. Two-centered orbitals 



Historically, the first use of two-centered orbitals for molecular calculations is attributable to Wallis and Hulbert's 
paper published in 1954. For a historical summary, see [33,E3|- New push and progress have been made by a 
school of French researchers fH fl^l fl^ flHl F^l since 1976. Their work has made the two-center orbital approach 
to diatomic modeling and computation an admirable success. 

The contributions by the French school are manifold, generalizing most of the aspects of one-centered orbitals to 
two-centered ones. To introduce its basic elements, let us utilize some of the semi-tutorial material from Scully, et al. 

^ . . . 

Recall the correlation diagram for H2 molecule shown in Fig. 1181 17]. A typical form of the trial wave function for 

H2, e.g., may be represented as 

tp = [ci'tijH2,lsa + C2'4'H2.2pa]f{ri2), (VI.52) 

where ^H2,isiy and ipH2,2pcr are chosen to be, respectively. 



ihi2M<y = V'ij+4^(l)V'i/+,i<,(2), 
i'H2ap<y = V'i/+,2p(l)V'//+,2p(2)- 

^From HVI.52|) and HVI.53|) , we see that the trial wave function HVI.52P 



(VI.53) 



(i) is uncorrelated if f{ri2) = 1- If /(''12) has explicit dependence on ri2, then (|V1.52p is correlated. A simple 
choice of f{r 12) was 

/(n2) = l + ^ri2, (VI.54) 

in Scully, et al. [T^ . but many such correlation functions / satisfying the interelectronic cusp condition as given 
in Subsection IvIB may be used also; 

(ii) has configuration interaction if C1C2 7^ in (|V1.52|) . 



It is found by [13 that 



58 



(a) for uncorrelated orbitals and without configuration interaction, i.e., f{r 12) = 1 and C2 = in (|VI.52p . the choice 
of 

i^H+Mij) = A^e-"^^^- [1 + i?2P2(M,)], J = 1, 2, (VI.55) 

m (|VI.53|l i. where P2 is a Legendre polynomial, gives the binding energy Eb = 0.132 a.u.= 3.59 cV for H2. 
Ifere, A/" is a normalization constant, and ai and B2 are two variational parameters. 

(b) for correlated orbitals but without configuration interaction, i.e., f{ri2) = 1 + ^ri2 and C2 = in ljVI.52|l . the 
choice of (IVI.55II yields the binding energy Eb = 0.1710 a.u.= 4.653 eV for H2. The choice f{ri2) = 1 + Kri2, 
where k is a variational parameter, slightly improves the answer and gives Eb = 0.1713 a.u.— 4.661 eV. 

(c) for correlated orbitals with configuration interaction, i.e., f{ri2) = 1 + ^?'i2 and C1C2 7^ in ljVI.52|l . the choice 
of HVI.55P along with 

^Ht .2P (^■) = (^^- ) + B^P^ )] (VI.56) 

in l|VI.53|l (ai, B2, a2 and B3 are variational parameters) yields the binding energy Eb — 0.1712 a.u.= 4.658 
eV. For the correlation factor /(r'12) with the variational parameter k we obtain Eb = 0.1721 a.u.= 4.682 eV. 

What is quite striking here is that if, instead of (|VI.55p . we choose a finer two-centered orbital 

V'i/+,i.(i)-AAe-"i^^[l + B2P2(M,) + S4P4(/^,)], J = 1,2, (VI.57) 

and perform calculations using three variational parameters ai, B2 and B4, then numerical results manifest that B2 
totally dominates -64, with a ratio |i34/i?2| ~ 2 x 10^^. This shows that the simple orbital (jVI.52p with (jVI.55|l is 
able to capture the essence of chemical bonding of H2 with very good accuracy. 

For a heteronuclear molecule such as HeH+, a simple adaptation of the above scheme works equally well. Never- 
theless, the authors have found that for large i?, more terms involving both A and pi variables should be included in 
(|VI.55p . such as 

^H+Mii)^^^~""'''[^ + B2P2{^iJ) + ■■■ + B2,nPn^{H)][l + A^L^{xJ) + ■■■ + AMx3% (x, =2ai(A,-l)), (VI.58) 

where Ln{x) are Laguerre polynomials, so that good accuracy can be maintained in the calculation of E. 

Another choice of simple two-centered wave functions, in the same spirit of this subsection, was given by Patil |57j . 
where he has generalized the one-centered Guillcmin-Zener type one-centered molecular orbitals (see Item (3)) by 
considering the gerade state 

= iV(l + bXfe-"^ cosh(a/i) (VI.59) 

and ungerade state 

= N{1 + bXfe-"^ sinh(a^) (VI.60) 

for (or any homonuclear) ionic orbitals. Asymptotic behavior of H^-like orbitals can be built in through /3, which 
plays the same role as /3 in the Jaffe solution ljIV.22|) . 

(1) Le Sech's simplification of integrals involving cross-terms of correlated wave functions 

Siebbeles and Le Sech developed an ingenious approach to the evaluation of energy by applying integration by 
parts (or, Green's Theorem) to avoid the quadratures of cross terms of the type 

V,($,$j) • V,/, i,j = l,2, (VI.61) 

where / is a correlation fmiction, and ^i^j is a product of two-centered orbitals. This is a fine feature of their 
approach. 

They choose wave functions of the form 

V'(ri,r2) = 0(1,2)17(1,2), (VI.62) 
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where f2(l,2) plays the role of the correlation function (but may be more general than the) / discussed in (|VI.52|I . 
while 



(l,2) = $(l)<i>(2) 



(VI.63) 



where $(i), i = 1,2, are the Hylleraas-type two-centered orbitals; cf. ljVI.77p below. Note here that each of <I>i and $2 
can have higher order molecular configurations such as 2s<Tg, SpCTu, 3dag and 4/cr,( (in the computation of He, e.g., in 

0)- 



Let the diatomic molecule's Hamiltonian be 



H=--{\7l + \7l) + V{l,2), 



where 



(VI.64) 



l^(l,2) = y(l,2) + — 
ri2 



with 



y(i,2) = - 



ria rib r2a 7'2b 



Let $„(i), i — 1,2, be eigenstates of the two center problem 



1 , = 1,2. 



Denote 0^(1,2) to be solutions of 



-2(v? + ^2) + ni,2) 



(/.,(l,2) = i?,,^,(l,2) 



Then 



1 



0,(1,2) ^ <J>„,(1)$„,(2) • ^[a(l)/?(2) - a(2)/?(l)] 



(VI.65) 



(VI.66) 



(VL67) 



(VI.68) 



(VL69) 



satisfies (|VL68|I with 

Sj-e«i+e„2- (VL70) 
Now consider a trial wave function for ljVI.64|) in the form 

0(l,2)f^(l,2) 

where 0(1, 2) is of the form IIVI.69|I while ri(l, 2) is intended to model the Coulombic repulsive effect from the l/ri2 
term and, therefore, plays a similar role as (but may be more general than) the correlation function jirn)- Without 
loss of generality, 0(1, 2) and il(l, 2) are assumed to be real. 
Denote the 6-dimensional Laplacian 



Consider the matrix element 



i/y EE (0,l]|--V2 + y(l,2)|0,i^) 



dridr2 



--0,riv2(0,r!) + 0,0,r!V(i,2) 



dridr2 



(r2>,vg0j- + 0,0ji7v^r2 + 20ii]V60j • Vef^) + 0^0J^^^F(l, 2) 



Ti 



(VI.71) 
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To treat 7i, write 

and note through an easy verification the foUowing: 



Vein"" 



(VI.72) 



But the LHS of l|VI.72|l is equal to twice of Ti. So we now substitute one-half of the RHS of l|VI.72p for Ti in (|VI.71|I . 
obtaining 



dridr^ 



(VI.73) 



where the divergence term V •[•••] in (|V1.72|I disappears after integration in the 6-dimensional space provided that 
(^i, 0j , V(ri^) and $1^ decay fast enough as |ri|, \r2\ oo. The RHS of (|VI.73|I is further simplified by utilizing 



and by combining T2 with terms therein, leading to 



dridr2 <{ i( 



Now, with (p satisfying ljVI.68(l we obtain 



dridr2 



-(t)i(j}j\VeQ\ 



(VI. 74) 



(VI.75) 



Comparing (jVI.75p with (jVI.71|l . we see that all the cross-derivative terms Vg'/'j • Vgil have been eliminated, and the 
only "burden of differentiation" has been placed only on iVgf^P in ljVI.75p . 

If we choose n = l + ^ru, then iVe^Jp = |VifJp + \\/2^\'^ = 1/2 and (|VI.75|) becomes 



dridr2 



(VI.76) 



But VL doesn't have to be chosen as O = 1 + |fi2- Le Sech has preferred a special form of Q. as given in (|33.(2)). 
Note that tg and e„ can be used as variational parameters, so can be d in 103.(2)). In a personal communication 
from Le Sech to Scully (9/6/2004), Le Sech pointed out that the trial wave functions in the form of (|V1.62() with 
ri(l, 2) chosen as (03.(2)) is particularly suitable for the diffusion Monte Carlo method 01 ■ 

Regarding the two-centered orbitals $ni(l) and ^"712 (2) in HVI.67P and HVI.69p . Aubert-Frecon and Le Sech 0, 
(4)] used a truncated summation from the Hylleraas series solution, cf. HIV.30|I : 



K 



N 



'^> = ^\,fi,cj,) = ^ /rFr(M,0)E^«^"'^^"'^[2p(A-i)r/2£-(2p(A-i)), 



(VI. 77) 



Tl = 



where p is used as a variational parameter, and /™ and c„ are coefficients which can be determined from the Killingbeck 
procedures [TSf. 



(2) Generalized correlated or uncorrelated two-centered wave functions 
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The two-centered orbitals we have been using in this section to build up the molecular orbitals, whether they 
be correlated, uncorrelated, with or without configuration interaction, have been heavily influenced by the classic 
solutions due to Hylleraas and JafFe, cf. pV.30|l and l|IV.26|) . These two famous solutions were derived during the 
1920s and 1930s with great ingenuity, their greatest advantages being that a 3-term recurrent relation of the series 
coefficients. The 3-term recurrences further lead to continued fractions which have enabled mathematicians to perform 
asymptotic analysis. That was during a time when no electronic computers were available and it was perhaps the 
only way to perform any theoretical analysis at all. Nowadays, fast computers are readily available so we don't have 
to rely overly on 3-term recurrence relations. Five-term recurrence relations are just as good and can be treated with 
relative ease also by the Killingbeck algorithm [t^, for example. 

Recall that for the single-electron, two-centered heteronuclear problem, separation of variables leads to 



[(A2-1)A']' + 

[(1 - ^,^)M']' + 



A2-1 



1 - /i2 

The appearances of the above two equations suggest the ansatz 



A = 0, (i?i = R{Za + Zb)/2) 

M = 0, (R2 = R{Za-Zt,)/2). 



(VI.78) 
(VI.79) 



A(A) = £ /fc,iPr(A), Mil,) = £ A,2Pr(A^), 



(VI.80) 



k=0 



k=0 



where P™(A) are the associated Legendre polynomials. Substituting (jVI.80|l into (|VI.78p and ljVI.79|) and equate 
coefficients of P™(A) and P™(/x) to zero, see Appendix L, we obtain 5-term recurrence relations 



Akjfk~2,j + Bkjfk-i,j + Ckjfkj + Dkjfk+i,j + Ekjfk+2j — 0, k — l,2;j = 0, 1, 2, 



(VI.81) 



where 



A-kl — Ak2 

Bkj : 



p2(fc — ui — l)[k — m) 
(2fc-3)(2A:-l) ' 

2RAk-m) 

' J = 1,2, 



2fc- 1 



Cui^Ck2 = -k{k + \)-A- 



p^{k — TO + 1)(A; + TO + 1) p^ik — m){k + m) 



Dkj - - 



2Rj{k + m + 1) 



Eki — Ek2 — 



2k + 3 

p^{k + m. + l){k + m + 2) 
(2fc + 3)(2/c + 5) 



(2A: + l)(2fc + 3) 
J = 1,2, 



(2fc+ l)(2/j- 1) 



(VI.82) 



The two-centered orbitals derived above differ from those of Hylleraas and Jaffe. The ansatz (IVI.80|I is not the only 
way to obtain two-centered orbitals. It is possible to derive several other solutions in different forms using different 
orthogonal polynomials. Numerical results indicate that these two-centered orbitals here give accuracy compatible to 
that corresponding to Hylleraas and Jaffe type solutions. 



(3) Numerical algorithm 



In contrast to the explicit analytic formulas for the quadratures of James-Coolidge wave functions given in IV II P 
and in the affiliated Appendices J and K, here the integrals will be computed using the Gaussian quadrature routines. 
(The analytic formulas given earlier in I VH P can then be compared with those obtained here as a good check.) For 
the two-electron problem, the number of coordinates are six, so that the energy calculation requires the 6-dimensional 
integration. Using the cylindrical symmetry, we can reduce two angular variables to one. 

The next simplification is the choice of the wave functions. If we only restrict the wave functions to be the 
James-Coolidge type, the 5-dimensional integration can be divided as one two-dimensional integration and two three- 
dimensional integration. However, if we include the exponential function of inter-electron distance, that is e~'^^'^'^ , this 
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5-dimensional integration should be evaluated as multidimensional integration, and the numerical accuracy should be 
carefully investigated. 

Here, we only restrict the wave functions to be the James-Coolidge type. The usual integration has the following 
form. 

Ai(Ai)A2(A2)Afi(^i)M2(/i2)$(0i - <^)2)r'l2dXld\2d^l^d^l2d(t)^d(t)2 (VI.83) 

Using the r^2-identity, cf. (|VI.31|) . the exponent of ri2 can be reduced into or -1. For n = 0, the whole integration 
is divided into only 5 one-dimensional integrations. However, for l/ri2, by the Neumann expansion the form of wave 
function includes sum of Legendre Polynomials with proper coefficients; cf. (|VI.39|I . Especially, the form of the A-part 
is 

P(A<)Q(A>) = I "iJ^^J^J^^J' \' I ^ (VI.84) 

This makes the form of integration over the A variable as 

/•OO /.OO 1*00 pXi pOC yOO 







(VI.85) 



The implementation of these numerical integration is automated by standard computer software. Let us describe 
the numerics for a homonuclear dimer H2 molecule and a heteronuclear dimer HeH"'' molecular ion. 

The approximation of wave function of two-electron system into a multiplication of one-electron system can be 
easily made in the natural orbit expansion 

vl/(ri, r2) = CfcXfc(ri)Xfc(r2) (VI.86) 

where Xfc(ri) is the wave function of one-electron two-center problem. 
For the ground state, it's well-known that ci ~ 1, or 

*(ri,r2) ~xi(ri)xi(r2). (VI.87) 

There are a few alternatives to approximate the xi(ri). One is using the exact solution of one-electron two-center 
problem, IsCg state. The Jaffe's form is 

Xi(ri) = e-"^(1 + Ar^.g„f^j ^/2,nP2™(M), (VI.88) 

cf. (II V. 2611 . 

And the other is Patil's wave function 

Xi(ri) = (1 + bXfe-"^^ cosh(a^); cf. ^VTE^. (VI.89) 
We may also include configuration interaction, such as 

^'(i"i,i"2) = cixi(ri)xi(r2) +C2X2(i"i)x2(r2) (VI.90) 

where xi is IsCg state and X2 is 2pCT„ state. 

The result is shown in Fig.^l Before the diagonalization (that is, CI), the asymptotic behavior of the ground state 
E{R) is monotonically increasing. By diagonalizing, the behavior at large R becomes almost flat, however, E{R) is 
slightly above the exact asymptotic value, -1 (htr). 

Next, we perform computations with correlation by using /(ri2) = 1-1- ^fi2. Improvement can be readily seen in 
Fig. 1201 which is much closer to the exact calculation done by Kolos 9]. 

In the calculations below fFigs. [H^ a nd I20|l . "exact" solutions of H^-solutions were used (whose coefficients are 
obtained through truncated Killingbeck [73 procedures). Instead, one can use (|VI.55|) . (|VI.56p . (|VI.57p by optimizing 
the coefficients a and B therein, or, for the Patil's wave function, by optimizing the coefficients a and f3 in ljVI.59p 
and ^VI.60|) . By doing so, we obtain the energy curve of the ground state as shown in Figs. 1211 and 1221 respectively. 
The energy curves are improved over a wide range of R values, and we obtain the binding energy of 0.171 a.u.= 4.65 
eV, close to the experimental value. However, E{R) in Fig. [2] overshoots the dissociation limit. 
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01 23456789 10 



R(a.u.) 

FIG. 19: Comparison of potential curves of the H2 molecule computed with no correlation factor using the exact solution of 
the one-electron wave function IVI.88II (small dots) and the Patil's wave function (IV1.89t (solid line). The inner (outer) curves 
are the result without (with) configuration interaction. Squares are the "exact" ground state E{R) from Ref. |^. Upper curves 
correspond to "excited states". 




01 23456789 10 



R(a.u.) 

FIG. 20: Same computations as with Fig. 1191 but with the correlation factor. 

VII. ALTERNATIVE APPROACHES 

A. Improvement of Hartree— Fock results using the Bohr model 

The Bohr model can also be applied to calculate the correlation energy for molecules and then improve the HF 
treatment. Fig. [23 shows the ground state potential curve for H2 molecule calculated in the Bohr-HF approximation. 
Such an approximation omits the electron repulsion term l/ri2 in finding the electron configuration from Eq. (jl.2|) . 
The difference between the Bohr and Bohr-HF potential curves yields the correlation energy plotted in the insert of 
Fig. E3I 
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Wavefunction 

'I'{1,2) = 1'(1)^'(2){1 + Kr^2) 
where 1 electron wavefunction is 

) = e-«^ (1 + (X) + ^(x)) (1 + P^(ii) + P^{^)) 
witli Variational Paramtersa, A,, A^, B^, B,, k 

1 2 2 4 



Minimum E = -1.17214 (a.u.) at R = 1.4(a.u.) 



-1.05 - 



-1.15 - 



R(a.u. 



FIG. 21: Ground state potential energy curve of the H2 obtained using the truncated exact wave functions of one-electron 
system with the correlation factor and several variational parameters. The curve yields the binding energy of Eb ~ —0.1721 
a.u.= 4.684 eV. Squares are the "exact" ground state E{R). 




FIG. 22: Ground state potential curve of H2 computed from the Patil's wave function with the correlation factor and variational 
parameters. The curve yields the binding energy of Eb = —0.1713 a.u.= 4.662 eV. Squares are the "exact" ground state E{R). 
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FIG. 23: Ground state E{R) for the H2 molecule in the Bohr and Bohr-HF models. Insert shows the correlation energy as a 
function of R. 



In Fig. I^we draw the ground state E{H) for the H2 molecule obtained with the Heitler-London trial function 
that has the form of the combination of the atomic orbitals [68l |: 

* = C {exp[-Q;(rQi + rfc2)] + exp[-a{rbi + ra2)]} , 

where a is a variational parameter. Addition of the correlation energy from Fig. 1231 improves the Heitler-London 
result and shifts E{R) close to the "exact" values. The improved potential curve yields the binding energy of 4.63 eV. 



B. Dimensional scaling 



Dimensional scaling offers promising new computational strategies for the study of few electron systems. This is 
exemplified by recent applications to atoms, as well as and H2 molecules D-scaling emulates a standard 

method of quantum chromodynamics |lOj |. by generalizing the Schrodinger equation to D dimensions and rescaling 
coordinates The D ^ 00 limit corresponds to infinitely heavy electrons and reduces to a classic electrostatic 

problem of finding an electron configuration that minimizes a known effective potential. 

We start from the Schrodinger equation = £'4' for a two particle wave function ^'(ri, r2). The H2 Hamiltonian 
in atomic units reads 

^ - -^V? - + V{pi,p2, Zi,Z2, 0, i?), 

where the Coulomb potential energy V is given by Eq. I|I.1|I . In a traditional dimensional scaling approach the 
Laplacian is treated in D-dimension and the wave function is transformed by incorporating the square root of the 
Jacobian via ^E" J~^/^$, where J = {pip2)^~^ {sin (f>)^^'^ in cylindrical coordinates, and <j) is the dihedral angle 
specifying relative azimuthal orientation of electrons about the molecular axis |7|| . D-scaling in spherical coordinates 
is discussed in Appendix On scaling the coordinates by and the energy by 1//^, with /=(£>— l)/2, the 
Schrodinger equation in the limit L* — > 00 leads to minimization of the effective potential [t^ 

E=l + -^+Vipi,p2,Zi,Z2,^,R). (VII.l) 

2 \Pi P2J sm^0 
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R, a.u. 



FIG. 24: Ground state energy E{R) of the H2 molecule in the Heitler-London method and the improved E{R) after the addition 
of the correlation energy. Dots are the "exact" result from |^. 



The obtained electron configuration is sometimes called the Lewis structure because it provides a rigorous version 
of the famiUar electron-dot formulas introduced by Lewis in 1916 77]. Fig. 1251 (upper curve) displays the Z? — s- 00 
potential curve E(R) that exhibits no binding and substantially deviates from the D — 3 "exact" energy (dots). In 
the limit D ^ 1 the dimensional scaling reduces the Hamiltonian to the delta function model [t^ 



A variational solution of the one dimensional wave equation [79j is pictured in Fig. 1251 flower curve). Also shown is 
an approximation to E{R) for D = 3 obtained by interpolating linearly in 1/Z? between the dimensional limits: 

EsiR) = '^E^iR) + ^Ei{R). (VIL3) 

The interpolated D = 3 curve exhibits binding which indicates the feasibility of extending dimensional interpolation 
to molecules. 

A proper choice of scaling can improve the method. Here we discuss a dimensional scaling transformation of 
the Schrodinger equation that yields the Bohr model of H2 in the limit of infinite dimensionality ,jy. For such 
a transformation, the large-D limit provides a link between the old (Bohr-Sommerfeld) and the new (Heisenberg- 
Schrodinger) quantum mechanics. The first-order correction in 1/D substantially improves the agreement with the 
exact ground state E{R). 

In the modified scaling the Laplacian depends on a continuous parameter D as follows 

1 d f jj_^d\ 1 92 9^ 



For Z) = 3, Eq. (IVII.4II reproduces the 3D Laplacian. On scaling the coordinates by the energy by 1/ P (recalling 
f — {D ~ l)/2) and transforming the electronic wave function by 

* = (piP2)"'''-''/'$, (VIL5) 

the Schrodinger equation is recast as 



{K + U + V)<P = E<^>, 



(VII.6) 
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R, a.u. 

FIG. 25: Ground state energy E{R) of the H2 molecule in _D = 00, D = 1 and three dimensional interpolation (solid lines). 
Dots show the "exact" D=3 energy from Q. 



where 

Q2 Q2 Q2 Q2 / 1 ;^ X Q2 



^ {D-iy \dpl ' dpi ' dzl ' dzl ' \pI ' pI) d(j)^ 



{ D-2){D-A) 

2{D-IY \pI ' pI 



f/= Vfl_fnil_^(^ + J_). (VII.7) 



In the D 00 limit the derivative terms in K are quenched. The corresponding energy Eao for any given internuclear 
distance R is then obtained simply as the minimum of an effective potential 



^ = o ( — + — ) +Vipi,p2,zi,z2,^,R). (YIU 



2 [pl + Pi 

This is identical in form to that for the Bohr model, Eq. Ijl.2l) . and we thus obtain for Eoo{R) the same solutions 
depicted in Fig. |2| This result for Eoa{R) differs in an interesting way from that obtained in a traditional study of the 
D —^ 00 limit for the H2 molecule 76). Here, in order to connect with the Bohr model, it is necessary to incorporate 
only the radial portion of the Jacobian in transforming the electronic wave function via Eq. HVII.5|) . The customary 
practice, which employs the full Jacobian, introduces a factor of 1/ (sin^ (j)) into the centrifugal potential, as seen from 
Eq. (jVII.l|) . The modified procedure yields directly a good zeroth-order approximation. 

The ground state E{R) can be substantially improved by use of a perturbation expansion in powers of l/D, 
developed by expanding the effective potential of Eq. IVILSI) in powers of the displacement from the minimum ; 
for He and H^ this has yielded highly accurate results Terms quadratic in the displacement describe harmonic 
oscillations about the minimum and give &\/ D correction to the energy. Anharmonic cubic and quartic terms give a 
contribution. For the He ground state energy (the R — limit for H2) a first-order approximation yields ^, Jl| 

^,(0) = ^ (1 - ^1 . (VII.., 



[D-iy \ D J 

where E^o = —3.062 a.u. is the Bohr model He result 0,113. For D ~ '3 Eq. I|VH.9|) improves the ground state 
energy of He to E{0) = -2.906 a.u., which differs by 0.07% only from the "exact" value of -2.9037 a.u. pl |. 
To evaluate the 1/D correction for arbitrary R, it is convenient to introduce new coordinates 

Zl = -^(Zl - 22), Z2 ^ --^[Zi+ Z2). 



68 



The effective potential of Eq. ljVII.8|) then has a minimum at pi = P2 = Po, (p = tt, Z2 = 0. Along the coordinates pi, 
P2, Z2 and (j) the potential has a single well structure no matter what the internuclear spacing R is. However, along 
the 5i direction at R= 1.2 the potential changes shape from a single to a double well; such symmetry breaking is a 
typical feature exhibited at large D '84'|. To avoid the inaccuracy of approximation by a single quadratic form one 
can solve the Schrodinger equation numerically along the zi direction for the exact potential and add contributions 
from the harmonic motion associated with the other coordinates pi, p2, Z2 and 4>. The result is shown in Fig. 1261 
0. The 1/-D correction improves E{R) and predicts the equilibrium separation to be i?e — 1.62 with binding energy 
Eb = 4.38 eV. 




FIG. 26: Ground state energy E of H2 molecule as a function of internuclear distance R calculated within the dimensional 
scaling approach (solid lines) and the "exact" energy (dots). 



VIII. CONCLUSIONS AND OUTLOOK 



Many major topics on diatomic molecules, and some other atoms and molecules in general, have been addressed 
in this article, giving it a very "locally diverse" or perhaps a somewhat disjoint look. But a simple, "global" picture 
may be viewed and understood in/from the following diagram: 
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old Bohr(-Sommerfeld) 
quasi-quantum-mechanical model: 
properties, new progress and refinement 
(Chapters 1 and 2) 



new Schrodinger(-Bom-Oppenheimer) 
fully-quantum-mechanical model: 
one- and two-centered orbitals, 
asymptotics, theoretical and numerical 
methods: Chapters 3, 4, 5 and 6 



"Unification" through D-scahng, 
and further (Bohr-Hartree-Fock) 
refinement: Chapter 7 



This "unification" is by no means an easy task. Nevertheless, it was a goal we somehow envisioned to achieve when 
this research was started and we hope we have made at least some success. 

At present, in the study of ultrafast laser applications to chemical physics and molecular chemistry problems, there 
is a pressing need to understand the quantum- dynamical behavior of molecules and the associated properties of excited 
states and their computations. Many challenges are lying ahead and awaiting elegant solutions. 



APPENDICES 

APPENDIX A: SEPARATION OF VARIABLES FOR THE H^-LIKE SCHRODINGER EQUATION. 

Let us consider (|IV.13|) . Here we show how to separate the variables through the use of the ellipsoidal (or, prolate 
spheroidal) coordinates (see Fig. [SJ 



R 



R 



2 V(A2-l)(l-/i2) cos,^, y=-^{X^-l){l-^^) sine 



z = — An. 
2 ^ 



(A.l) 



Note that the coordinates A, n and are orthogonal, and we have the first fundamental form 



dx^ + dy"^ + dz'^ = hi dX^ + hi dfi'^ + h 



where 



hi 



K = 



dx 
dX 
dx 

dx 



dy_ 
dX 

dy_ 

dy 



dz 
dX 
dz 
dfi 

dz 



R^ l~fi^ 
R^ X^-l 
R^ 



Thus 



1 



h\hf,hr,: 



d_ 

dX 



h^h^ d 
hx dX 



h\h^ d 



d 



(A2 - 1) 



d_ 
dX 



d_ 

dfi 



d^ 

(1-M^) 



d ( hxhf, d 



dfi 



(A2-l)(l-^2) 



(A.2) 
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Note that through the coordinate transformation (jA.ip . we have 

equivalently, 



A = 
M = 



R ' 

ra - n 



R 



Also, we have A>1,— l</i<l. 
Write 

* = A(A)M(/i)$(0). 
4>((/>) must be periodic with period 27r. Therefore 



$(0) 



imtj> 



m = 0,±l,±2,. 



Substitute (|A.2|I . (|A.4p and (|A.5(I into (|IV.13(I . and then divide by e 



1 



d_ 

2i?2(A2-^2) 
2 ^6 



Af - 



,2^ 5 



A 



(A^ - fi^)m^ 
(A2-l)(l-/i2) 



MA 



-MA 



2 ^" MA+^AfA-SA/A. 



RX + fi i? A - ■ i? 

Further muhiplying every term by — -^(A^ — /i^), we obtain 



d_ 
dX 



d 



A- 



A2-/i2 



(A2-l)(l-^2) 



RZhiX -fi) + RZaiX + /i) - ( - ^ ) (A' - m') 



RZaZf) R^E 



MA = 0. 



Set 



P 



-(-R^E + RZaZb) > 0. 



(A.3) 



(A.4) 
(A.5) 



(A.6) 



(A.7) 



(A., 



We have > here due to the fact that we are mainly interested in the electronic states that are bound states, i.e., 
not ionized. 

Let the constant of separation of variables be A. Then from ljA.7|l and (|A.8|I we obtain ljIV.161) and l|IV.17p in 
Subsection CvlB. 



APPENDIX B: THE ASYMPTOTIC EXPANSION OF A(A) FOR LARGE A 



Here, we provide a quick proof of Eq. (|IV.22p . From pV.19(l . we have 



A"iX)~-p^A{X) 



2A 



A^-l 



A'(A) 



A + 2RiX 
A2 - 1 



A2-1 



(A2 - ly 



A(A) 



(B.l) 



We now find the Laurent expansions of the coefficient functions on the right-hand side of IjB.ip as follows 



2A 



A2-1 



2 

A ^ A2 



2 

A3 





A4 



2 

A5 



A + 2RiX \ , ^ / 1 1 

^^-^(A + 2i?,A) 1 + - + - 



2i?i A 2Ri 
A2 +~ 



A_ 

A4 



71 



2 



A2-1 ; A A2 A3 A4 

m2 2rn2 



± 



(A2-l)2 A A2 A3 A4 AS A6 
^From these expansions above, we now use the ansatz 

A(A)=aoe-^'^A'^ (l + £l + £| + . . . ) 

by substituting it into ljB.l|l and equating aU the coefficients of A^" to zero for n = 0, 1, 2, . . . . We easily obtain 

P = 1, ci = 



p 2Ri-2pl3-p 

All the other coefficients c„ can be determined in a straightforward way. Note that the asymptotic expansion ljIV.22|) 
also gives 

n 

A(A) - aie-P^X^ ^2^^ 0{e-P^X^^"+^^), for A > 1. (B.2) 

APPENDIX C: THE ASYMPTOTIC EXPANSION OF A(A) AS A ^ 1 

Here we provide a proof of Eq. (|IV.23p . Multiply HIV.18|) by (A - 1)/(A + 1) and rewrite it as 

■(A + 2i?iA-p2A2 



2A 

= (A-1)^A"(A) + ^(A-1)A'(A) + 



-(A-1)- 



A(A) 



A + 1 ^ ' (A + 1)2 

2 

« (A - 1)2A"(A) + (A - l)A'(A) - ^A(A), for A « 1. (C.l) 

A differential equation (such as (|C.l|l ) set in the form 

(x ~ lfy"ix) + (x - l)q{x)y'{x) + r{x)y{x) = 0, (C.2) 

near x = 1, where q{x) and r{x) are analytic functions at x = 1, is said to have a regular singular point at x = 1. The 
solution's behavior near x — \ hinges largely on the roots v of the indicial equation 

v{v - I) + q{l)v + r{l) = Q (C.3) 

because the solution y{x) of HC.2|I is expressible as 

oo oo 

y{x) = h{x - ly^ ^^(^ " + ^2(2; - 1)"' ^'^■(■^ - 1)'' = ^0 = 1)' 

k=0 k=0 

where 1^1 and 1^2 are the two roots of the indicial equation l)(13p . under the assumptions that 

vi > 1^2, ^1 " ^2 is not a positive integer. (C.4) 

However, if (|C.4|I is violated, then there are two possibilities and two different forms of solutions arise: 
(a) v\ = V2. Then 

y{x) = biyi{x) + 621/2(2;), (C.5) 

where 

00 

^/l(x) = (x-l)''l^Cfc(x-l)^ (co = l) (C.6) 

k=0 
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and 

oo 

2/2(:r) = (x-l)''^^dfc(2:-l)'= + [ln(a;-l)]yi(a;); (di = 1). (C.7) 
fc=i 

Solution 7/2 in (|C.7p should be discarded because it becomes unbounded at x = 1. 

(b) i/i — 1/2 = a positive integer. Then case (a) holds except with the modification that 

oo 
fc=0 

(do = 1, c is a fixed constant but may be 0). (C.8) 

Applying the above and ljC.3|) to equation IjC.ip : 

2 

977 ^ 

(A - lfh"{X) + (A - l)A'(A) - — A(A) « 0, (C.9) 

we obtain the indicial equation 

o 

777 

+ — = 0, (C.IO) 

with roots 

Vi = \m\/2, v-i — — |to|/2, (m can be either a positive or a negative integer). (C.ll) 

Thus either 

v\ = V2 (when m = 0) 



v\ ^ vi = \m\ ^ a positive integer, where m 7^ 0. 

Again, we see that solution 2/2 in HC.8|) must be discarded because it becomes unbounded a.t x = 1. Thus, from HC.6|I 
and HC.9|I . we have 

00 

A(A)«(A-l)^^c,(A-l)'=. (C.12) 

APPENDIX D: EXPANSIONS OF SOLUTION NEAR X^l AND A » 1: TRIAL WAVE FUNCTIONS OF 

JAMES AND COOLIDGE 



In the pioneering work of James and Coolidge :20j , the two-centered trial wave functions for H2 are chosen to be 
(ground state) 

^(Ai,A2,Mi,M2,p) = ^e-"(^^+^^' J2 C„„,fep(ArAViM2+A'/A™/4Ai^2)/, (D-1) 



2tt 



and 



(excited state) 



V'(Ai, A2, Ml, M2, P) = e-"(^^+^^) J2 CrnnMK^2^^{t^2 - A? A^^'/.J/i^^)/, (D.2) 
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Where p = ri2 is the distance between the two electrons of H2, the coefficients Cmnjkp are computed from the 
minimization of the energy, plus possible orthogonality conditions. 

We examine the asymptotics of the trial solutions (|D.1|I or l)D.2(l based on the discussions in Subsection IIV Bl We 
see that the exponent e^""^^ in e~"('^i+'^2) would correspond to the factor e^^^ in (|IV.2ip or ljIV.26|) . This is excellent 
as it reflects the exponential decay in the radial variable of the (first) electron. However, the other non-constant 
polynomial terms of the form 



^1 ^2 



\n2 \ m2 k j 
A2 MlA*2' 



with \mi\ + \ni\ > 1, 



^=1,2, 



(D.3) 



possess polynomial growth rates in either Ai or A2, which are at odds with the asymptotics in IIIV.22p since for large 

"^Aa' or Af A"" in XPlk besides the 
(One might argue that the polynomial growth A7'^A2^ or A"^A™^ would be 
-p^2 xhis can be true, however, only by increasing p and thus it causes the 



A (which may be either Ai or A2), there should not be any polynomial growth A^^Aj^ or A"^A 



exponential decay factor e 



-pAi 



-pX2 

-pXi 



killed by the exponential decay e 
loss of accuracy.) One might still argue that the typical term in the series (jIV.26|l for m = (ground state) satisfies 



(A + 1)^ 



A- 1 
A + 1 



= A'^ 



((T- l){a-k - 1) /I 

r2 V A 



for A > 1. 



This means that either mi or ni in l|D.3|l should never exceed a. Can't we, at least, use terms A™^A2^ 
some restrictions such as 



or A^'^A^^ with 



< mi, m < a; 



(T = — - 1, for TO I , ^=1,2? 
P 



The most important behavior of ^E* happens near A = 1. For Awl, the typical term in ljIV.26|) satisfies 



(A2-1)V(A + 1)° 
2^+--'=(A-l)^ + 



A- 1 
A + 1 



|TO| 

2 



a — k 



A - 1 



\_m[ 

2 



a-k-1] 



1 • 2 



= (A-1)^[2- 



(A - 1)'^ + ©((A - 1)'^'+^)], for A « 1, A > 1. 



(D.4) 



When TO = 0, the case of the ground state, the above expansion in terms of powers of A is consistent with the terms 
involving powers of Ai or A2 in l|D.3|l since those powers in l)D.3|l can always be re-expanded in terms of powers of 
Ai — 1 or A2 — 1. However, when to = 1 (or to = ± odd integer, for that matter) for the excited state, the function in 

(|D.3|) is no longer consistent with those in ljD.41) as far as the A- variable is concerned because the factor (A — 1)t^ in 
(|D.4(I is unaccountable for those in (|D.3(I . Indeed, we have indicated in (|IV.23(I through asymptotic analysis that the 

factor (A — is inherent in the solution and, thus, must be properly taken into account. 

The discussion in this section points out some weaknesses in the choice of basis functions IjD.ip or ljD.2p based on 
the asymptotic arguments for A 3> 1 and A ~ 1. Such weaknesses may have contributed to the fact that many terms 
are required in (ID.1() in order to calculate or to do the variational analysis of the energy E accurately for H2 by James 
and Coolidge ,2Qj . 

Our conclusion for this subsection is: because of the vastly different asymptotic behaviors of A(A) for Awl and 
A 3> 1, the best strategy for numerical computation is to use two different representations for A(A), one for A « 1 and 
another for A ^ 1 and match them, say, at a medium-size value such as A = 5 or 10. This, however, will invoke more 
computational work and is beyond the interest of the authors for the time being. 



APPENDIX E: THE MANY-CENTERED, ONE ELECTRON PROBLEM IN MOMENTUM SPACE 

Here we derive the eigenvalue equation ljIV.43p . We start from the model equation 




■0(r) = Eijjir) 



ix,y,z) £ 



(E.l) 
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Recall the Fourier transform 



VI/ (p) 



(27r)3/2 



Note that the Fourier transform of the potential term 

1 



(27r)3/2 



\r-R,\ 



ilj{r)dr. 



7dr 



(E.2) 



(E.3) 



is a divergent integral in the classical sense. However, the modern mathematical theory of the "regularization of 
divergent integrals" |23, Chap. 3] makes (|E.3p well defined: 



= (2^)-3/2 J e-^ir' +Ri)-pl_dr' {r'=r- Rj) 
= (2^)-3/2e-««. P lim / e-^*- P- dr' 



(27r)-3/2e-^«rP lim / / p ^ - r'^ dr' sin 6 d6d<j) 



47r 



a^O p2 ^ Q,2 



2- (27r)-i/2g-«fi,-P/p2^ 



(E.4) 



Applying the Fourier transform to (|E.ip by utilizing (jE.4|) and other well known properties such as convolution, we 
obtain 



N 



-iflj-(p-p') 

———^.(p')dp' = 0. 



Define 



then we obtain the integral equation 



N 



\P-P 



(E.5) 



(E.6) 



(E.7) 



Now, we project the 3-dimensional momentum vector p onto the surface of the unit sphere, S^, of the 4-dimensional 
space, the 1-1 correspondence ^ p through 



= '^PoPxip'^ + pI)~^ = sin X sin 6* cos 0, 
^2 = 2poPy{p^ + pI)-^ = sinxsin6'sin(/), 
6 = '^PaPzip'^ +Po)"^ = sin X cos 0, 
^4 = {pI - P^){P^ + Ply^ = cosx 
O<X<7r,O<e'<27r,O<0<7r, 



} ^eM^^ = 1^1 = 1, 



while keeping in mind that 

Px = psin0cos(/), py — p sin 9 sin (j), pz=p cos 6. 
Then for ^ ^ p, ^' ^ p' , and k be the angle between ^ and we have 

COSK = ^ 

1^ - = ^2 ^ ^'2 _ 2^ . ^' = 2 - 2C0SK = 4sin2(K/2). 



(E.9) 



(E.IO) 
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Also, it is straightforward to verify that 



Hence 



1 



\p-p\ 



{p'' + pIW^ + pI) 



pI 

{p^+pUp'^+Pl)s\u'{nl2y 



Let <Kl be the infinitesimal surface area of the 4-dimcnsional hypersphere 53. Then from ijE 

dTL = SIT? X sin ^ dxdOdcj) = — sir? x dxd{cos 

while IjE. 9(1 gives the standard 

dp = P^ sin 6 dpd9d(f) = —p^ dpd{cos 

^From I|E.8P a. we further obtain 

dp _ p'^ + pI 
dx 2po 



The equations ((E.8(l and l|E.12|l - ((E.14|l now give 



dp = — 



p^ dp 



sin^ X dx 



dn = 



? I 9 \ 3 

P +Po 
2po 



dn. 



We can now use (|E.11|I to eliminate the denominator inside the Fourier integral of l|E.7|) . Moreover, define 
where g ^3 is the point ^ ^ p. Then the integral equation ((E.7|l simplifies to 



Po ^ J 8tt^ sin^(K/2) 



tp{n')dn' 



(E.ll) 

(E.12) 
(E.13) 

(E.14) 
(E.15) 
(E.16) 

(E.17) 



'S3 



At this point, we need to introduce the hyperspherical harmonics on O3, which constitute an orthonormal basis for 
square summable functions on 53 and are given by 



Y^emin) = {-iYCi{x)yirn{0,4>), 

n = 1, 2, . . . , ^ = 0, 1, 2, . . . , TO , + 1, . . . , 0, 1, . . . , £, 
cf. jEHl-lEni) for YimiO,(t)), 



(E.18) 



where 



cUx) 



2n{n~-e~ 1)1 



TT{n + £y. 



-|l/2 f 


■ d ' 


(sin^x) j 


d(cosx)_ 



C„_i(cosx) > , 



(E.19) 



are the associated Gegenbauer functions, and the C^'s are the Gegenbauer functions, with the generating function 



^/i^Q(m), \h\<l. 



(E.20) 



For^,f e 



1 



1 



'|2 



ei-2cos«;(e/c) + (e/0' 
1 1 



ife7e< 1, 



(E.21) 



.c'M -2 cos + 



7T2 ife/?'<l- 



Denote ^> = max(^,^'), ^< = min(^,^'). Then ljE.20|) and (|E.2ip give the Neumann expansion 



l\2 



^ CO 

7T X! 



^> n=0 



C„(cos k) 



n=l ^> 



which, in the Umit as — > 1, yields 



4sin2(K/2) 

OO 

= ^C„_i(cosk). 



(by m^)) 



n=l 



But by the addition theorem of angles, we have 



27r 



rt ^ — ^ 



so from (IE. 2211 we obtain 



n> 1,^== 0,1,2,..., m= -£+!,..., 0, 



1 °° 1 



and, thus, the kernel of the integral equation ljE.171) can be written as 



where t = (n£m) runs triple summation indices according to (|E.24|) . We now use the orthonormal basis 
(jE.lSp to make a re-expansion 

e'^^'-^Ytin) - J2 snRj)Yrin), J = 1, 2, . . . , iV, 

r 

where r = (ji'i'm') runs triple summation indices similarly to t. Then 



and IE. 2511 becomes 



SiiRj) = J e'^^-PYtin)Yrin)dn 



Z,e-^^r(p-p')/lSn'sm'i>,/2)]^Z,Y.[SnR,)r^Sf{R,)YAmr^^ 



The orthonormal expansion of (^(il) in (|E.17p is denoted as 



Substituting ljE.27|l and ljE.28|) into (jE.lTp and equating coefhcients, we obtain 



PQ 
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This is an eigenvalue problem 

Pc = poc (E.29) 

where P is an infinite matrix with entries 

p}, ^j2^^J2iSrim*}/riR,) 

j T 

where c is an infinite-dimensional vector with entries Cf . 

The value of po will yield the energy E from ljE.6|) . One can obtain the wave function V-'(^) from applying the 
inverse Fourier transform to ^(p) through (|E.16|I and (|E.8(I . 



APPENDIX F: DERIVATION OF THE CUSP CONDITIONS 



Here we derive the cusp conditions at the singularities of Eq. ljV.3|l . Since the idea is the same for each of the five 
sets of singularities, we will only treat the most complicated case, ri2 — 0. 

Use the center-of-mass coordinate system, see Appendix O we can transform (|V.2|I into the form 



2Za 



2Zb 



2M ^ 2fi ■'■1^ |2-R„+ri2| \2Rh+ri2\ 
'2Za 2Zfc q\q2 ZaZ\, 



|2i?a-ri2| |2i26-ri2| 



r\2 



R 



(F.l) 



where S is the center of mass coordinate of two electrons, see Appendix [U] for the details of the notation. 
We now define the spherical means of a function. Given a point Tq € M^, the spherical means of a function u{r) at 
ro on the sphere with radius p is defined to be 



^(ro) = ^JJ u(ro + pv)du 



(F.2) 



where Sp is the sphere with radius p centered at rg, here p = 1; diu ^ sin 9 dOdcf), where lo represents all the angular 
variables; v = (vi, 1^2, 1^3), + V2 + 1% = 1] v is the unit outward pointing normal vector on Si. It is easy to see that 
if u{r) is continuous in a neighborhood of Tq, then the spherical means just converge to the pointwise value: 

lim Uav,p{rn) = w('"o)- 

p-tO 

We now integrate over a small 3-dimensional ball Bp^ with radius pQ centered at ri2 = 0, for any S E M^, S 7^ 0: 



1 1 07 



V'(ri2,5) 



2Zb 



\2Rt+ri2 
'7192 



■V'(n2,5) 



2Za 



\2Ra~ri2 



2Zt 



\2Rb-ri2 



■^(ri2,5) 



ri2 



-'ip(r 



12 



5) + ^V'(n2,5)|rfri2=0. 



Note that as po 0, the integrals of all the terms above vanish (by their continuity at ri2 = 0), except possibly those 
of 



-;lv?2^(ri2,5), ^V^(ri2,5). 
2fi ri2 



Thus, we need only consider 



lim 

po^O 



-^V?2V'(r-i2,5) + — VXri2,5) 
ri2 



dri2 = 0. 



(F.3) 
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Apply the Gauss Divergence Theorem to the first term of the integral to get 



ri2 



ip(ri2,S)dri2 



(F.4) 



But dSpg = Pq duj; thus 



dip 



dSon 



dtp 
dr 



12 



. 2 2 11 dip{ri2,S) 2 



dr 



12 



-duj = 47rpQ 



dms) 



dr 



12 



(F.5) 



where Si — Spg\pg=i. By using spherical coordinates, dri2 = dujdri2, we have 



glg2 
''12 



ip(ri2, S)dr 12 = J J J '^^tp{r 12, S)di 



r\2dri2 



where e(ri2,S) — > as ri2 0. This follows from the fact that ip{ri2,S) is continuous at ri2 — for any S. 



(Continuing from the above) = (71(72 



rPo 

ipiO,S) / ri2dri2 +£(po,S')47r 
"'0 



(7i(72VX0,5)47r • ^ + 2nple{po,S). 



iFrom ijO)) . l|R4jl . (jF3|l and (ESjl, we now get 

"9V'(0,5)1 



1 2 



9ri2 



au,po 



-^p^^i(72V(0,S) + |e(po,5). 



Dividing all the terms above by p§ and let pQ — > 0, we have £(po: iS) ^ and, therefore. 



lim 

Po^O 



a'/'(o,5) 



dri2 



= pqiq2ip{0, S), for all S € 



(F.6) 



(F.7) 



(This is equation (2.11) in Patil, Tang and Toennies ^3|.) Similarly, one can derive the electron-nucleus cusp condition 
at, e.g., ria = 0, for IjV.lj) to be 



lim 

Po^O 



di>{ri,r2) 



dria 



-miZaip{ri,r2)\ri^=o, for aU r2. 



Theorem F.l. Assume that mi — m2 and (71 = (72 = —1 in l|V.l(l . Let {S,ri2) denote the CM and relative 
coordinates and UJ12 denote the angular variables of the vector ri2. Let tp be a nontrivial solution of (jV.ip such that 
its local Taylor expansion near ri2 = is of the form 



^{ri,r2) = Co(S) + Ci{S)ri2 + ©(r?^), for r^ small, 



(F.9) 



where Cq{S) and Ci{S) are independent of oj 12 while the remainder satisfies 0(r^2) — C'3('*i2, 'S')?'^2 /'^'^ some bounded 
function C3 which depends on ri2 and S. Then Ci{S) — ^Co{S) for all S and, consequently, 

^(ri,r2) = Co{S) (1 + ^ru) + ©(r?^). (F.IO) 
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Proof. We substitute the RHS of (jF.9|) into (|F.7p . There is no need to take the spherical mean on the left of ljF.7|l 
anymore as the dominant term on the RHS of Hi^'.10|l does not depend on the angular variables of ri2. We therefore 
obtain 

Ci(5) = fiqiq2Co{S) = ^Co(5), as ^ = mim2 ^1 = = 
2 mi + m2 2 



Hence 



CoiS) 



l + iri2 



which is llFJl). □ 



It looks as though the condition (|F.9|I is somewhat contrived. However, useful application can be seen shortly, in 
Theorem El 

Similarly, we can obtain the cusp conditions at ria = 0, rif, = 0, r2a ~ and r2b = 0, as given in the following. 



Theorem F.2. Let tp be either a nontrivial solution of (jV.ip or a trial wave function for (|V.1|I . Let ria and LOia 
denote, respectively, the radial and angular variables of the vector ria- Assume that for ria sufficiently small, tp 
satisfies the Taylor expansion 

^(ri,r2) = Co(r2) + Ci(r2)ri„+0(r2j, r^a small, (F.ll) 

for some functions Cq and Ci depending on r2 only, where the remainder satisfies 0{r\^) = C'i{ri,r2)r\^ for some 
bounded function depending onri andr2- Then 

'0(^"l-^'2) = C'o(r2)(l - m^iZaria ) + 0(r2j, ria small, (F.12) 

where 0(r?„) = C3(ri,r2)r?<j. 

Proof. The kinetic energy term ^ 2mr^i' ^^er the translation 

R 

Zi zi + — , ci. IjCip and Fig. [^for notation, 
becomes — 2mr^ri„ which is centered at {xi,yi, zi) — (0,0 — -j). Now apply the Hamiltonian 



2mi ria 



1 ^2 -^6 Zb ^ qiq2 ^ ZaZh 

2to2 ^ rib r2a r2b ri2 R 



(F.13) 



to l|F.lip . We need only focus our attention locally near ria = and ignore the terms in the bracket of (|F.13|I as they 
have no effect on the singularity at ria = 0. We obtain 

(H - E)i, = --}-^J-[rl ■ Ciir2)] - ^[Co(r2) + Ci(r2)riJ 



2mirf^dria'''^ ria 
higher order terms in ria 
1 2Ci(r2) Za-Coir2 



— , higher order terms in ria- 

2mi ria ria 

To eliminate the singula^rity cit vi^ — above, it is necessary that 

^ + ZaCo{r2) = 0. 

mi 

The above gives (|i^'.12(l . as desired. □ 

Example F.3. The trial wave function ^p{ri,r2) — (f>{i'ia)4'{i'2b), where (j>(r) = e^"^ is a one-centered orbital, satisfies 
condition (IFTTIi of Theorem Ol □ 



One can easily apply Theorem lF.2l to other singularities at rib = 0, r2a — and r2b = 0. 
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2 1 2 2 



Z Z Z Z I 
ria rib 7'2b ri2 



1 



z^ 

R ' 



Theorem F.4. Let 4>{r) he a sujficiently smooth function. Let the Hamiltonian represents a homonuclear case: 

H 

Then the product function 



1 + -^n. 



satisfies the interelectronic cusp condition as given in (jF.9|) . 

Proof. We represent the variables ri and r2 in terms of the CM coordinates 



r2 ==5'- i(ri -ra) 



where S = \{ri+ r2). 

Then for ri2 small, by Taylor's expansion, 



^•2 



5 + -n2 



5--n2 



V^iS) ■ ri2 
V(/)(5) • ri2 



1 1 
2!4 
1 1 
2!4 



Therefore 



ri)(?!)(r2) ( 1 + iri2 



1 1 
2! 4 
1 1 



</.(5) + ^V<^(5) • ri2 + ^-r?2 • D^^iS) ■ r^^ 



q^{S) - ^V0(5) • 7-12 + ^ • D^q^iS) ■ ri2 ± 



1 + 2-12 



^^(iS) ( 1 + 2 ''12 ) + quadratic or higher order terms involving ri2. 



Hence condition IF. 91) is satisfied 



□ 



APPENDIX G: CENTER OF MASS COORDINATES FOR THE KINETIC ENERGY - j^V? - a^Vi 



The coordinates for electron 1 and 2 are, respectively, 

ri ^ {xi,yi,zi), r2 = (a;2, 2/2,22)- 

The kinetic energy operator is 



where 



(G.l) 



1 -v2 1 V72 



2mi ^ 2to2 



^2 _ ^ SP_ 
^ dx^- dyj dz'j 



J = 1,2. 



Define the CM (center-of-mass) coordinate S: 



TOiri + m2r2 

mi + 7712 
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I 

-RJ2 R/2 

FIG. 27: Various vectors are defined in this diagram. 



and (relative coordinate) 

ri2 =ri -r2. 

Here we derive the kinetic energy term in coordinates S, ri2. For any difFerentiable scalar function /(ri,r2), we have 

df = Vi/ • dn + V2/ • dr2 = Vs/ • dS + V12/ • dri2, 
where Vs and V12 are the gradient operators with respect to the variables of S and ri2. Then 



But 



Hence from HG.2p and (|G.3p . 



[Vl/,V2/] 



dri 
dr2 



[Vs/,Vi2/] 



dS 

dri2 



dS 




dri2 





mi m2 
mi+m2 mi+m2 

1 -1 





dri 




dr2 



[Vl/, V2/] = [V5/, V12/] 



mi+m2 mi+m2 

1 -1 



(G.2) 



(G.3) 



1 

2m 1 





1 

2m2 



Vl 

V2 



[Vs, V12] 
[Vs, V12] 



mi+m2 mi+m2 
1 -1 





2m 1 



2m2 , 



mi+m2 
mi+m2 





"Vs" 




Vi2_ 







2(mi+m2) 

_if^ + ^) 

2 \^ mi m2 y 



Vs 
V12 



H 



1 



1 / 1 



2(mi+?7l2) 2 \TOl 7712 

1 9 1 9 
Vo V?n, 

2M ^ 2^ 



where M = mi + 7772 is the total mass of electrons and fi 



mi m2 
mi+m2 



is the reduced mass. 
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APPENDIX H: VERIFICATIONS OF THE CUSP CONDITIONS FOR TWO-CENTERED ORBITALS IN 

PROLATE SPHEROIDAL COORDINATES 

In the work of Patil (see Eq. (2.15) in he indicated that for the ground state ip (i.e., with azimuth quantum 

number m = 0) for the molecular ion with Hamiltonian 

the "coelescense" condition at r;, = can be expressed as 

hi, \ 

= ~Zb, (H.2) 

rb=0 

the angle 9 is introduced in Fig. |2H1 He then indicates that the above is "essentially the same" as Kato's cusp 
condition. 



(x,y,z) 



/I di^ 




1 d^j 




\4' dn 


H 


"0 drb 






-R/2 R/2 

FIG. 28: A local spherical coordinate system (ra,S, <^). 

Actually, at least two ways are viable, which are going to be described below through some concrete examples. 
The first way takes a limiting approach in a similar spirit as Path's (|H.2p . The second way utilizes an explicit 
representation. 

Example H.l. For the Hamiltonian (|H.1|I . motivated by the Jaffe solution PV.26|I . let us consider a trial wave function 
for the ground state: 

V'(A,A*) = e-"^[l + B2F2(A*)]. (H.3) 



We want to consider the cusp condition in terms of the undetermined coefficients a and B2 in ljH.3|) . 
Recall from (IA.3|I that 

(i) ?'a — > is equivalent to A ^ ^ —1; (H.4) 

(ii) r;, — + is equivalent to A ^ 1, /i ^ 1. (H.5) 

The singularities in the Laplacian in prolate spheroidal coordinates, according to (|A.3p . are discerned to be 
contained in 



— ji^ X + fi A — /.i 



Thus, we deduce that for (jHTTJ, after substituting HH.3|) into H, the following: 

(i) 

HiP r-. Fi(A,^), for « 0, where A > -1, (H.6) 

A + ^ 

where we have dropped terms not containing the singularity (A + ^)~^ and collected the dominant terms corresponding 
to singularity (A + /i)^^ in 

Fi{X,pi) = -J^J^^—^ [^2aXe-"\l + B^P^iii)) + (A^ - l)a'e-''^{l + B^P^^^jl)) + e-"^3i32(l - 3^')] 

9 7 

-^e--'[l + B,P,{^,)] 
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(ii) 



where, similarly, 



Hi; ^ -^Gi{X,ii), for rf, w 0, where A > 0,/L(« 1, (H.7) 
A ~ /i 



For ljH.6p and ljH.7p to stay bounded, it is necessary that Fi{l, —1) = and 01(1, 1) = 0. 

1 17 

-1) = -j^ h2ae-"(l + B2) - 6i?2e-"] - -^e-"(l + ^2) = 0, 



I.e., 



I + B2 

Gi(l, 1) = - [-2ae-"(l + B2) - 6^26-"] - -^^-"(1 + B2) = 0, 

i.e., 

These are the cusp conditions at ra = and rb = 0. 

However, we need to remark that 1) = and Gi(l,l) = are only necessary conditions for the desired 

boundedness because the limits in HH.6() and HH.7() do not exist in general as A and /U may be related in infinitely 
many different ways to yield totally different limits as A — > 1 and fi ±1. Thus, the above estimation approach has 
an ad hoc nature. Nevertheless, (|H.8|I and (jH.Op do provide correct answers as to be cross-validified with ljH.16|l and 
(pLT7)l in Example HOI 

In inspecting l|H.8|l and l|H.9p . we see that they are consistent when and only when Za = Zb, i.e., the homonuclear 
case. Therefore, (|H.3|I would not be a good choice of a trial wave function in the hcteronuclear case. For the 



heteronuclear case, taking hints from the exact solutions (|IV.30p and (|IV.35p . we choose the trial wave function 
Here, a may be different to (3 even though the exact solution says a = j3. We have 



A + /X 



where the singular terms are collected in 



F2{\ m) = - ^2(;f_^,) {-2aAe-"^e'5^[l + B,P^{y) + i?2P2(M)) 
+ (A2 - l)a2e-"Ag,3M(i ^ B^P^{^i) + B2P2{ti)] 

+ e^"^e''^[(/3' - - + B^P^{^i) + B2P2{^^)) + (2/3 - 2/3/i' - 2m)(Bi + 3^2^) 

+ (1 - ^l^)iB2)] - — ^ie-"^e'^^[l + B,P,{^i) + B2P2{fi)], 
H 

with the terms corresponding to the dominant singularity in 

F2(l, -1) = - ■^e-'^'^ [-2a{l - Bi + B2) + (2/3)(l - Bi + B2) + 2(Bi - 3S2)) 
"^e— '^(l-i?i+i?2), 
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a — (3 = RZa 



(Bi - 3^2 ) 
I-B1+B2 



Similarly, the behavior near = gives 



where 



A — /i 



G2(A,/z), 



G2(A,m) 



-{ - 2aAe-"^e''^[l + BiPiifi) + B2P2{^i)] 



i?2(A + ^) 



-QA^/3Mr//32 _ ^2,2 



ip' - /3'^i' - 2/3/i)(l + BiPi(p) + S2P2(m)] + (2/3 - - 2f^)iBi + SSa/x) 
2Zb 



(1 - /i2)3B2]} - -^e-"^e'^^[l + + i?2P2(A*)]; 



G2(l, 1) = - (-2a(l + Si + S2) - (2/3)(l + Si + ^2) - 2(Si + 3B2)) 

-^e-"+^(l + i?i+S2), 
it 



I.e., 



a + (3 = RZb - 



(Bi + 3S2) 
1 + Si + B2 



So, the cusp conditions give 

Bi = 

B2^- 



3{RZa + 2/3 - RZb) 



2(3 + 4a - 2i?Za - 2RZb - RZaU - RZaP + R^ZaZb - aRZb + (3RZb + - P^) ' 
(2^2 + 2a- 2RZaa - 2aRZb - RZa + 2PRZb - 2/32 _^ 2R^ZaZb - RZb ~ 2RZa(3) 



2(3 + 4a - 2RZa - 2RZb - i^Z^a - RZafi + R^ZaZb - ai?^^ + l3RZb + a^ _ /32) 
If, in addition, we set a = /3, then 

3{RZa + 2a- RZb) 



Bi = 
B2 = 



2{2RZaa - R^ZaZb + 2i?Za - 4a + 2i?Z6 - 3) ' 

[-2R^ZaZb + RZa + 4i?Z„a + - 2a) 
" 2{2RZaa - R^ZaZb + 2RZa - 4a + 2i?Zb - 3) ' 



□ 



Example H.2 (Verification of the cusp conditions via explicit calculations). Consider the same Hamiltonian as given 
in Mil. 



We translate the origin from (0,0,0) to (0, 0, -~R/2) and set up spherical coordinates (ra, 0, (f) as shown in Fig. |2H1 
^From the cosine law, cf. Fig. [2H1 

rb = {rl + i?2 - 2Rra cos 9)^''^ 



„ , , 2 cos 1 2 



1/2 



i? 



COS0 2\ 

^-^ra + 0{rl) 



(H.IO) 



As we are exploring the singularity behavior of Hil) near Va = 0, we need only concentrate on the dominant terms of 
Ta and, thus, we drop 0{r^) in IjH.lOp and approximate by 
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Using ET3|) . we therefore obtain 



A= — (l-cos6i) + l, 
R 

^=!^(l + sine')-l. 



(H.ll) 



The transformation l|H.ll|l will greatly facilitate our calculations. 
Consider the trial wave function l)H.3(l again here. We have 



1 - 



Bo 



+ 3B2H 



dra 



(H.12) 



^From l|H.ll|l . we have 



dX 


1 


dra 


" R 


dfi 


1 


dra 


" R 


JH.12| 


as 



A - 1 

ra 
Ai + 1 



(H.13) 



dra [ 

+3B2{{p + 1) - I) 



A- 1 

M+ 1 



f g-"(A-l) 

ra 



-[a{l + B2){X-l) + 3B2{pi + l)] 



-a • ^2 • 2(A - + 1) + 3B2(Ai + 1)^ 



-B2-a{X-l)ifi + lf 



(H.14) 



If we set ra — p and then substitute l)H.14(l into the left-hand side (LHS) of IjF.Sjl . the two rightmost brackets on 
the RHS of ljH.14p will become 0{ra) = 0{p) and C(r^) = 0{p'^), and then vanish after p ^ 0. So we only need to 
consider the spherical mean of over a sphere with radius p which is given by 



47rp^ J p 
5„ 



~Q (1 — COS 6) 



^(1 + E2) • 4 (1 - cos 0) - 3^2 (1 + cos 61)1 p2 duj 
R Hi 



27T fTT 



4vr ./o Jo 



COS 9 



a(l + B2)--+3B2 
H 



-a(l + i?2)-|+3i?2--^ 



cos U > sm 



1 _Q,_ii£ 



1 



^)[a(l + B2) + 3S: 



2j ■ 



R 



-^(e'=-e-'=) + i(e'= + e-'^) 



[-a(l + B2) + 3B2] 



1 

R 



(H.15) 



with fc = ap/ R. Letting p ~+ 0, we have fc — > 0, and 



(RHS) of JILT5ll = -e-" • + B2) + 3B2] ■ - 

R 



(RHS) of (|F.8|I (modified for the case of molecular ion) 
-Z,e-"(1 + B2), 
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a(l + B2) + W2 = RZa{l + B2), 
which is exactly (jH.8|l . Similarly, at = 0, we can obtain 

a(l + B2) + 3B2 = RZbil + B2), 

which is exactly (|H.9|I . 



(H.16) 
(H.17) 



APPENDIX I: INTEGRALS WITH THE HEITLER-LONDON WAVE FUNCTIONS 



The integrals involved in Eq. IjVI.SI) can be separated into eight types of elementary integrals: 
(i) I J a(l)(-V?)a(l) dx(^^\j b{2){-^l)h{2) dy ^\ j a{2){~^l)a{2) dy 

R3 R3 R3 

= \j &(1)(-V?)6(1) dx^=^-, 

R3 



(ii) \ J a(l)(-V?)5(l) dx(^=l J a(2)(-V?)6(2)d. 

R3 R3 

(iii) J a^{l) (j-:^^ dx = -a; 



'x = 1 + w w' 



(iv) S = / a(l)5(l) dx^e-"^ [l + w + 



w 



(v) aJ= / a\l) ( -— \ dx = a 



(vi) aK = J a{l)b{l)(^--^^ da; = -ae-'"(l + u>); 



(vii) aJ' ^ II a^{l)b^{2) [ — ] dxdy = a 



-2w 



1 11 3 1 . 

\ 1 — w H — 

w 8 4 6 



(viii) aK' = / a(l)6(l)a(2)6(2)— dxdy = -al - e-^- (^_^ + + + iu;3 

ri2 5 I \ 8 4 3 



6 



+S'^E^{-Aw) - 2SS'E^i^2w)] }, 



where 



w = aR; 

7 = Euler's constant = / dt - dt = 0.57722 • • • ; 

Jo t Ji t 

00 

f e"* 

Ei{x) = integral logarithm ~ — (P.V.) / — ^ dt (for x > 0), 



here P.V. means "principal value" of a singular integral. 
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APPENDIX J: DERIVATIONS RELATED TO THE LAPLACIAN FOR SUBSECTION IVIlD 



Evaluation of the Laplacian in the prolate spheroidal coordinates for the general form of the wave function which 
includes electron-electron correlations explicitly is a demanding task. Here we provide the details of the calculations. 
The general form of each term in the wave function is 



$,(1,2)- 2^<f,(l)<I>,(2)P,(ri2) 



with 



(J.l) 



and 



$,(2) =F,(A2)Af,(/X2). 



It can readily be seen that the Laplacian operates only on <&s(l)- 
According to l|A.6|l : 



V2$,(l) 



4 



(A?-M?) 



(A? - 1)(1 - g^2 + Q^^i^i - 1)— + —(1 - Ml) 



dXi d^i 



dm 



{Xi-^if)F,Gs d^Ps 



(Af-l)(l-A.?) 



d 2 ^ MF^Ps 



We first single out part of the second term in the square bracket above for further evaluation to obtain 



d 



d{FsPs) _ d 



dXi 



{K - 1) ^p. 



'dXi 



dPs 
'dXi 



ox<''-'^ ox. 



Similarly, part of the third term inside the square bracket of ljj.4p becomes 



^'dX.^^'-^'dX, +^^^i"^^9Ai 9Ai 



F ^ (X^ 11^^^ 



d , 2 



d{GsPs, 



d 



9^1 



dPs 



—{l-^Jii)lPs^ + Gs^ 



d ^dGsdPs 



(J.2) 
(J.3) 



(J.4) 



(J.5) 



(J.6) 



We can rewrite the Laplacian to take into account the separable functional dependence of the three parts of the 
wave function, that is. 



1 



$s(l) 



V2$,(l) 



VfF, , V?G, , V?P, , 2ViP, -ViP, ^ 2ViG, -ViP 



P, 



G, 



P,P, 



+ 



GsPs 



(J.7) 



since 



Vi$,(l) 



2(A?- 1)1/2 ^ 2(1-^2)1/2 g^^(^) 



P(A2-Ai?)i/2 dXi 
2 



R{xj~^4y/^ a/ii 

9$.(1) 



i?(A2- 1)1/2(1 -^2)1/2 --PI' 

where, 6^1,6/^1, and e^j are the unit vectors pointing to the respective directions. 



(J.8) 
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Now, wc use the actual functional forms of various parts of the wave function to complete the evaluation of the 
expectation value of the Laplacian. Setting i^s(Ai) = e~"^^X^% we obtain 



d 



R^{\l-iil)d\x 
4 d 

^2 Ti2\ 



R^{\l-Hi) 
+m,e-"^i ((m, + 1)A7« - (m, - 1)A™=-^) 
4F, 



m,ae-"^i(A7»+i - A7»-i)} 



i?2(Af-Mf) 
Similarly, setting Gs(/Ui) = we obtain 



a (Ai - 1) - 2a (m^ + l)Ai - nis— + mg (tos + 1) - (m^ - l)--2 



Ai 



'A? 



4 a 

AGs 



(1 - 



(J.9) 



Since the effect of the Laplacian on a given function is coordinate-free, we can consider the evaluation of Vf Ps{i'i2), 
through the spherical coordinates for a general function f{r 12) of ri2 to obtain 



V?/(ri2)=V^^J(ri2) 



1 d 2 #(''12) 



ri2 dri2 dri2 



Thus, with Ps(ri2) = ri2, 



V?P, = — {rlAl - M2' + 2r-i2M2 = ^^P. 
ri2 ri2 



(J.IO) 



(J.ll) 



The other terms involved are 



It' 



.idr\2 



dXi ^'12 9Ai 
RHP, 



8r?2 
RHP, 



2Ai - 2A2/U1M2 



Al - A2/Ui;U2 - 



M 



Xl-l 
MAi 



COS(01 - (j)2)'^Xl 



cos(^i - ^2) 



P2;p^ Al 

4r?2 (A?-l) 
P^^P, Al 
8r?2 (A?-l) 



(Af - 1) - A2/Ltl/i2 



(Af-1) 
Al 



— MCOS(01 — 02 ) 



r-2 + \2 _ ^ 2 _ 2 _ 2 , 2A2M1M2 
-^ri2 + Al - A2 - Ml - /i2 H ^^^^ — 



(J.12) 



and 



dPs RHP, 



RHP, 



2/^1 — 2A1A2/U2 



j--^ _^2-) ^os{(j)i - (f>2){-2ni) 



^'^12 



Ml - A1A2M2 + _ ^2^ COS(0i - (p2) 



RHP, Ml 

■ 4rf2 (1-M?) 
P^^Ps Ml 

■ 8r?2 (l-M?) 



-(1 - M?) + AiA2M2-^^-^ - Mcos((^i - 02) 
Ml 

4 2 ,2 %2 , 2 2 , 2A1A2M2 
^^12 - Al - A2 + Ml - M2 + 

-K^ Ml 



(J.13) 
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Upon putting the above together, the Laplacian operation part of the Hamiltonian takes the form 
1 _ / VfF. , VlGs , VlPs , 2ViF, • ViPs , 2ViG, • ViP, 



'i>.(l) 



P. 



R^{Xl-nl) 
4 



G^Ps 



a\\i - 1) - 2a (m, + l)Ai - m,— + mj (m, + 1) - (m, - 1) 



+ 



1 



MI 



(A?-M?) 
1 



(-aAi + TOs) — 



Ai 



4 2 , x2 i2 2 2 , 2A2AilM2 

^^12 + '^1 ^ -^2 ^ A^l ^ A*2 " 



A? 



+ 1) ) + 



Ai 



-^ri2 - Ai - A2 + Ml - /"2 



2A1A2M2 



Ml 



(J.14) 



As before we can construct the inter-term expectation value integral for the Laplacian using the above relation. 
Introducing the function 



and defining 
we obtain 



X^im.n.j.k-.l) = Z''{m,n + 2,j,k;l) - Z^im^nJ^k + 2;i), (J. 15) 

{■m,nj,k;i) = {mr^rirjr^kr'jr) + {ms,ns,js,ks]is), (J. 16) 

Ai Aa'Ml M2 ''12 I I Ai °A2'Mi M2 ''12/ 

: [a^X''{m + 2, 71, j, k- 1) + {-a^ + (m, - j,)(TO, + + 1 + 4) } X^im, n, j, fc; £) 

- 2a(7is + l)X''(m + l,n,j,k]i) + 2ansX''{m- l,n,j,k;i) 

+ m,(m, - l)X"{m - 2, n, j, fc; £) + - l)X''{m, n, j - 2, fc; ^) 

+4(4 + 1) {X^im, n + 2, j, k;£~2)- X^im, n, j, k + 2; £ - 2)} ~ 4aA^''(m + 1, n, j, k; £)] 



- gsaX^im + 3, n, j, /fc; ^ - 2) + 4aX''(TO + l,n + 2, k; £ - 2) 

+ isaX^im + l,n,j + 2,k;e-2)+ isaX^im + 1, n, j, fc + 2; ^ - 2) 

- 24aX''(TO, n + 1, j + 1, fc + 1; ^ - 2) 

+ 4(m, + js)X''{m + 2, n, j, k;£-2)- 4(™. - js)X''{m, n + 2, j, k,i-2) 

- 4(m, + js)X''{m, n, j + 2, A:; f - 2) - 4(m, - js)X^m, n, j,k + 2-J- 2) 

- 2lsmsX''{m. - l,n+ l,j + l,fc + 1;^- 2) - 2is3sX'' {m + l,n + I, j - 1, + 1; ^ - 2) . 

Thus we have furnished complete details of the electronic kinetic energy calculations. 

APPENDIX K: RECURSION RELATIONS AND THEIR DERIVATIONS FOR SUBSECTION IVllD 

In this appendix we provide simple proofs of the recursion relations which are needed in the analytical calculations. 

1. A{m\a) 



A{m;a) = j A™e-"^dA 



(K.l) 



m 

a ./I 



= — [e " + TO^(m - 1; a)]. 
a 



When TO = 0, 

A(0;a) = / e-"^dA 

Ji « 

The recurrence relation can be used in succession to give 

A TO;a — > 7 TT — ■ 

2. F{m;a) 

The definition is 

/oo 

To prove the recurrence relation for F{m; a), 

POO 

mFim - 1; a) - (m - 2)F(m - 3; a) = / (mA""^ - (to - 2)A™-3)g-aAg^^;s^)^;^ 



(A™-A'"-")e-"^go(A)|Q +ay dA (A™ - A™-^)e-"^Qo(A) 

/OO 
dA A"-2g-"A^ 

So, the recurrence relation is 

F(to; a) = -F(to - 2; a) + - [to_F(to - 1; a) - (to - 2)F(to - 3; a) - A(to - 2; a)] 

a 

The initial conditions for F(m; a) are already given in IIVI.44|) and l|VI.45|) . 

3. S(m,n;a) 



S{m, n; a) is defined by 



The recurrence relation is 



with 



S(m, n;a) — — [mS{m — 1, n; a) + A{m + n; 2a)] , 
a 



S(0,n;a) = -A{n-2a). 
a 



By using S{m, n; a), the following integrals can be represented in terms of A{m; a) and S{m, n; a): 



/oo /"Ai /"OO /"Ai 

dXiJ dA2ArA5e-"(^i+^^) = y dXiXTe-"^' J rfAaA^'e 



/ A" 



Ai 

n 

1 



+ ^ / dAsAr'e^"^^ 



^(!!i±Z!iM + £ZA(TO;a) + ^ TdAi /'^ dA.Ar Ar^e-"(^^+^^) 
a a a .h .h 



= > — Aim + n — v,2a)- H > —A[m\a) 

i/=0 ^ ' v=Q ^ ' 

= > — A{m + s;2a)— -{ > — --A(m;a) 

s=0 ^=0 ^ ' 

— — S'(n, m; a) + A(ri; a)A(TO; a). 
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Furthermore, 



/OO />C30 
dXi J rfA2A7*A^e-"(^^+^^) = A{m; a)A{n; a) 

/ dXi / dA2A^*A^e-"(^i+^=) + / dXi dAsA^^A: 



n -q(Ai+A2) 



So 



dAi / dA2A'i"A^e-«(^i+^=) + / dXi dX2y^X^e-"^^^+^^^ 
II Ji Ji Ji 

= — S{n, to; a) + A{n; a)A{m; a) — S{m, n; a) + A{m; a)A[n; a) 



S(m, n; a) + S(n, m; a) — A{m; a)A{n; a). 
4. T{m,n;a) 



The definition is 

Tim, n; a) ^ ^ ^F{n + v- 2a). 

The recurrence relation and T(0, n; a) are 



T(to, n; a) = — [mT{m — 1, n; a) + F{m + n; 2a)] 



a 



and 



T(0,n;a) = -F(n;2a) 
a 

5. i/o(n*»i;Q!) 

By definition, 

/OO />oo 
dAiy dA2A7'A^e-"(^i+^=)Qo(A>) 

dXi dAsAr A^e-«(^i+^=)Qo(Ai) + dXi dX2X^Xfe-''^^'+^'^Qo{Xi). 
The first term on the right-hand side above yields 

dXi [ ' dA2A7'A^e-"(^i+^^)go(Ai) 
1 Ji 



(K.IO) 
(K.ll) 

(K.12) 

(K.13) 
(K.14) 



(K.15) 



y°°dAiA7e-«^iQo(Ai) j 

i: 



dA2A^e-"^= 



dAiAre-"^igo(Ai) 1 -^e 



+ - 
1 ail 



dA2A^-'e 



n—l—a\2 



F(m + n; 2a) e " ^ , , n , , , 

— ^ ■ H F{m; a) + -Ho{m, n-l;a) 

a a a 

F{m + n; 2a) ^ e"" 

—5^ H F{m; a) 

a a 

nf F(m + n-l;2a) er'" , ^ n-l , 

+ ■ + F(m;a) H Ho(m,n - 2; a) 

a \ a a a 



1 \ - 1 , n! e " \ - 

= > — F(m + n- v]2a)-, -7 H > - 

a a" ^ (n- i^y. a ^ a 



(n — v)\ 

1 a^ ^ .n\ e~" v-^ 
-> — F(m + s;2a)— + > — 

s=0 



L/=0 



'{n-v)\ 



F{m; a) 



nl 



s\ a ^—ia'^(n — v)\ 



F{m; a) 



= — T{n,m;a) + A{n;a)F{m;a), 



(K.16) 



while the second term in ljK.15|) is the same as the first term if we interchange m and n. Therefore, 
Ho^m, n] a) = —T(n, m; a) + A(n] a)F{m] a) — T{m, n; a) + A[m] a)F{n; a). 

6. Hi(m,n;a) 

By definition, 



oo />oo 



Hi{m,n-a)= dXi dA2A5"A^e-"(^i+^^)Pi(A<)Qi(A>) 



Ai 



dXi / dA2A5"A^+ie~"(-^i+^^'Qi(Ai) + (same as left, with m n). 



The first term in l|K.18(l is evaluated as 

/oo /• Ai 

dXi J dA2A™A'2'+ie-"(^i+^^'Qi(Ai) 

/OO /'•^i poo p^i 

dXi j dX2X'^+^y^+^e-"^^'+^^-'>Qo{Xi)- J dXi j dA2A™A'2'+ie-"(^i+^^) 

dA2A™+iA'^'+ie-"(^i+^^)Qo(Ai) - 5(m, n + 1; a). 

By combining the two terms in ljK.18(l . we obtain 

Hi{'m, n; a) = HQ{m + 1, n + 1; a) — S{n, m + 1; a) — 3(171, n + 1; a). 

7. Hr(m,n;a) 

The recurrence relations for the Legendre polynomials are 



(r + l)Pr+l = (2r + l)xPr - TPr-l, 

(r + l)Qr+i = (2t + l)xQr - rQr-i. 

We then have 

Hr{m,n;a)= [ / A™A^P^(A<)Q^(A>)e-"(^i+^=) dAidAa 



AfA^' [(2t - l)A<P,_i - (r - 1)P,_2] 
[(2r ~ l)A>g,_i - (r - 1)Q,_2] e-"^^'+^'^ dX,dX2; 
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T'^Hri'm, n; a) 

= (2r - lfHr-i{m + 1, 71 + 1; q) + (t - lfHr-2{m, n; a) 

- (2t - l)(r - 1) J J \Tm\<Pr-iQr-2 + A>Pr-2Qr-i)e-"(^^+^^^ dAidAa (K.24) 
= (2t - lfHr-i{m + 1, n + 1; a) + (r - lfHr-2{m, n; a) 

~{2t-1) j j XTy^{X<{{2T - 3)X<Pr-2 - {r - 2)Pr-3)Qr-2 

+ XyPr-2{{2T - i)XyQr-2 - (t - 2)Q,_3))e""'^^+^'^ rfAidAa 
= (2t - l)2iJ^_i(TO + 1, n + 1; a) + (r - lfH-r-2{m, n; a) 

- (2t - l)(2r -i) I I A™A^(A'< + A|)P,_2Qr-2e-"(^^+^^) dAidAa 



+ (2t - l)(r - 2) / / A5"A^(A<P,_3Qr-2 + A>P,_2Qr-3)e""^^^+^^' dAidAa (K.25) 



(2t - lfHr^i{m + 1, n + 1; a) + (r - l)^iJ^_2(™, ^^; a) 
- (2t- l)(2r-3)(ff^_2(™ + 2,n;Q!) + iJ^_2(™, " + 2; a)) 



+ (2t - 1) y y A™A'2'(A<P,_3((2t - 5)A>Q,_3 - (r - S)Qr-i) 

+ A>((2r - 5)A<P.-3 - (r - 3)P,-4)Qr-3)e-"(^^+^^) dAidA2 
(2t - lyHr-iim + 1, n + 1; a) + (r - lfHr-2{m, n; a) 
- (2r - 1)(2t - 3)(i7r-2(rn + 2, n; q) + Hr-2{m, n + 2; a)) 

+ 2(2t - l)(2r - 5) / / A™A'2'A< A>Pr_3Qr-3e-"'^'+^'^ dAidA2 



(2t - l)(r - 3) / / A^ A^(A<P,_3Qr-4 + A>P,_4Qr-3)e-"(^^+^^' dAidA2. (K.26) 



As shown in the equations l|K.26p - (|K.26p above, the same patterns are repeated until r is reduced to 0. So, let's 
consider the last term. If r is even, the last term is 



-(2r-l) J y A™A'^(A<PiQo + A>Pogi)e-"(^^+^^)dAidA2 

= -{2t-1)J y"A™A'2'(A2<Qo + A^Qo - A>)e-"(^i+^^)dAidA2 
= - (2t- l){HQ{m + 2,n;a) + Hoim,n + 2;a)) 



dXi J 



dX2 A5"+U^e-"(^i+^=) 



+ (2t-1) / dXi dA2 A™A^+ie-"(^i+^=^) 

- (2t- l)[ifo(m + 2,7i; a) + iJo(TO, " + 2; a) 

- S{m + l,n;a)- S{n+l,m;a)] (K.27) 
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If T is odd, the last term then is 



(2r - 1) y y A^(A<Pogi + A>PiQo)e-"(^^+^=' dAidAs 

(2t - 1) / / ArA^(A<A>Oo - A< + A>A<Qo)e-"(^^+^=) dAidAa 



2(2t - l)iIo(m + 1, n + 1; a) 



g-a(Ai+A2) 



/CX) /■ Ai 

rfAi y dAa A7^A^+^ 

-(2r-l) / dAi / dA2Ai"+iA^e-"(^i+^=) 



= {2t -l)[2HQ{m+l,n+l-,a) - S{m,n+l-,a)- S{n,m+l-,a)\ (K.28) 

8. H^r\m,n;a) 

The recurrence relations for the associated Lcgcndrc polynomial arc 

{x" - l)^'^P!;+\x) = (r - z/)xP;(a;) - (t + //)P;_i(.t), (K.29) 
- lY'^Q^+\x) = (r - r.)a;g;:(x) - (r + ^/)g;;_i(x); (K.30) 

Hl^\m,n;a)=T^ J J X^X^iX^Pr - Pr-i){X>Qr - Qr-i)e-"^^'+^'UXidX2 (K.31) 

= T^[Hr{m + l,n + l;a) + Hr-i{m,n;a)] 

J J XTXq{X<PrQr-i + X>Pr-iQr)e-''^^'+^'^ rfAidAs 

= T^[Hr{m + l,n+l;a) + Hr-i{m,n;a)] 

- -Tz -T- [(2t + lfHr{m + 1, n + 1; a) + T^Hr-i{m, n; a) - (r + lfHr+i{m, n; a)] 

t(t + 1)2 t2(t + 1) 

= — —Hr+i{m,n; a) - t{t + l)Hr{m + l,n+ 1; a) H — —Hr-i{m,n;a). 

2r + 1 2r + 1 



So 



2r + 1 

i?^^)(m,n;a) = (r + l)i?x+i(m, n; a) - {2t + l)Hr{m + l,n + 1; a) + THr-i{m,n; a) (K.32) 



r(T + l) 

9. Hi^\m,n;a) 

Similarly to the preceding paragraph, 

= 11 ArA^[(Af - l)(Ai - l)Y/\{r - 1)X^P^ - (r + 1)P^_,) 

X ((r - 1)A>Q^ - (r + l)Q^_i)e-«(^i+^=) rfAidAs 
= [(r - l)2ijW(m + 1, n + 1; a) + (r + ^. 

- (r^ - 1) / / ArA^[(A? - l)(Ai - 1)]1/2(A<f;q1_i + A>F;_iQi)e-"(^^+^=^) dAidAs, (K.33) 
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and 



-ff|+i(m,7i;a) 



X ((2t + l)XyQl - (r + l)Qi_i)e-"(^i+^^) dAidAa 
= -^[(2t + (m + 1, n + 1; a) + (r + q,)] 



Therefore, 



(2t + ^. ^ ^2^^ _ l)H)^l^{m, n; a) - (r + 2)(r - 1)(2t + + 1, n + 1; a) 

+ (t + l)2(r + 2)iJ^L\(m,n;a). 



(K.34) 



(K.35) 



APPENDIX L: DERIVATIONS FOR THE 5-TERM RECURRENCE RELATIONS dVI.Sll l 



First, we cast both equations (|VI.78p and (|VI.79|) into the form 



0, J = 1,2, 



where 



Set 



fori = l, x = X, (f> = A{X), Ri = R{Za + Zb)/2; 
forj = 2, x^fi, <l> = M{fi), R2 = R{Za- Zb)/2. 



fc=0 



as in (jVLSOp and substitute HL.2|) into (|L.l 



oo 



dx dx 



+E^ 

fc=0 



-A-2RjX+p^x^ - 



k=0 

oo 



— (1 - x^)- - A 

dx dx 1 ~ x^ 



l-a;2 

Pk{^) + E A-[-2i?,x +p2x2]P,"(x) = 0, 



fc=0 



^ AhMfc + 1) - ^]Pr(^) - E hR,xP]:{x) + ^ hv^x^P^{x) = 0, 
^ h\~Kk + 1) - A]Pr(x) - ^ hR, ^fc-yKfc + "^)^ri(^) + (fc - ™ + 

i.-n fc=o 



fc=0 

oo 



fe=0 



(2/j + l)(2fc + 3) 



(fc — m + l)(fc + m + 1) (fc — m)(fc + m) 



^'^ ^"^^^ (2fc-l)(2fc + l) ^fe-2W^-0- 



(2A:+ l)(2A: + 3) (2/s - 1)(2A: + 1) 

We an now shift indices to convert Pfe"2(^)' -^fe" i(^) ^fe-2(^) Pk^i^)- We obtain 

— m — l)(fc — m) 



E 



fc=0 



rA;-2- 



2Pj(fc - m) 
(2fc-3)(2fc-l) 2fc- 1 



/fe-i 



fkCk 



„ 2P,(fc + m+l) p2(;;,_^„j^;L)(;;,_^j„_^2)' 

^, , „ 1- Jfc+2 zr; „s.^. J- = U. 



2A; + 3 



(2fc + 3)(2/c + 5) 



(L.l) 



(L.2) 
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The terms inside the parentheses above are exactly the 5-term recurrence relations (jVI.81|l . 



APPENDIX M: DIMENSIONAL SCALING IN SPHERICAL COORDINATES 

For description of diatomic molecules cylindrical coordinates provide a natural way of making a dimensional (D-) 
scaling transformation. Here we show how to do the D-scaling transformation in spherical coordinates, which is useful 
for description of atoms. Let us first consider the Laplacian in the Z?-dimensional hyperspherical coordinates 

xi = rcos^i sin6'2 sin^s • • • sin0£)_i, 

X2 — r sin^i sin02 sin03 • • • sin0£)„i, 

X3 = r 008 6*2 sin 03 sin 04 ••• sin 

X4 — r 008 6*3 sin 04 sin 05 ••• sin 0£)_i, 



Xj = r cosOj^ismOj smtljj^i ■ ■ ■ smtlD^i, 

xd-1 = rco8 0£i_2 sin0£i_i, 
xd = rco8 0£)_i, 

O<0i<27r, O<0j <7r, for j = 2, 3, • • • , L» - 1, (M.l) 

where _D is a positive integer and _D > 3. Define 

D~l 

h=l[h„ (M.2) 



where 



3=0 



D / ^ \ 2 
oxi 



Then the scaling factors are 

ho = 1, 

hi = rsin02 8in03 • • • sin0£)_i, 

/i2 = rsin03 8in04 • • • sin0£)_i, 

hk = r 8in0fc+i 8in0fc+2 • • • sin0z)_i, 

hD-2 = rsmOo-i, 

h-D-i = r, 

h = r'°-i8in02sin2 03sin^04---sin'="^0fe---sin'°"2 0£,_i. (M.4) 
The £)-dimen8ional Laplacian now becomes 



' D 



1 f 1 d . k-ln A_ 

r2 sin^ 9k+i sin2 0,+2 • • • sin^ Od-i I sin'^-^ 0,- 39^ ' 90, 



8in^-^0Z3_i— } (M.5) 
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Define the generalized orbital angular momentum operators by 

92 



^1 " der 



Li = - 



■ sm f 2 



sin 02 902 " dQi sin^ ^2 ' 



sm t/fe 



sm Wfc 



sin^ ( 



Then we have 



where 



Vl=KD-i{r) 



7^2 



KD-i(r) 



1 ^ ^.i^-i 5 



Let us consider Schrodinger equation for a particle moving in D-dimensions in a central potential V{r) 

To eliminate the angular dependence we separate the variables by writing 
Near the origin r = 0, 



with 



Vlr^Y{nD-i) = + D-2)- Cy-^Y{nD-i) = 

Ll_^Y{nD_^) = CY{nD_^). 



The effective Hamiltonian is given by 

Hd = -^Kn-i{r) + ^-^^ + V{r) 

With the following transformation 

=^-(S-l)/2$^ 

the corresponding equation for $ d reads 

1 52 A(A + 1) 



2 9r2 2r2 



V{t) 



$r> = Ed^d, 



where 



K = l + -{D-i). 

Equation IjM.lSp is the Schrodinger equation in D-dimensions for the function <f>£). 



(M.6) 



(M.7) 



(M., 



(M.9) 
(M.IO) 
(M.U) 
(M.12) 

(M.13) 
(M.14) 

(M.15) 
(M.16) 
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As an example, consider the Schodinger equation for the H-atom in D-dimensions 

l_o 



In the scaled variables 



3^ 

1 E 
2%' 



with tq — D{D — 1) / A and i?o = 2/(Z) — 1)^, the Schrodinger equation reads 

2 



D-l\ Z 



(M.17) 

(M.18) 
(M.19) 

(M.20) 



Now, let us write the Laplacian in spherical coordinates and transform the wave function 5" according to (|M.14p . 
We obtain 



1 / 3 



{D-l){D-i)\ i{D-l) Z 
2b V, 



* = E,^ 



2\DJ \dr^^ 4r2 

In the limit 13 — s- oo Eq. (jM.21|l reduces to a simple algebraic problem of minimization the expression 

9 5Z 



E, = 



8r2 2', 



(M.21) 



(M.22) 



which yields = 3/2Z and Es — —Z'^/2. This value coincides with the ground state energy of the hydrogen atom in 
3 dimensions. 
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